Step3_DE_analysis
Sam Pring, Sarah Lower, Brian Vestal
2025-03-12
Load packages
#installation if necessary
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#if (!requireNamespace("devtools", quietly = TRUE)) {
# install.packages("devtools")
#}
#BiocManager::install("rhdf5")
#BiocManager::install("tximport")
#BiocManager::install("DESeq2")
#BiocManager::install("tximeta")
#install.packages("tidyselect")
#loading in the data
#library(rhdf5) # allows for the handling of hdf5 file formats
#if (!require("tidyverse")) install.packages("tidyverse");
#library(tidyverse)
#library(tximport) # package for getting Salmon results into R
#library(datapasta) # allows for copy and pasting
#library(DESeq2)
library(SummarizedExperiment)
library(readr)
#library(tximeta)
#filtering and normalization
#BiocManager::install("vsn")
#library(edgeR)
#library(matrixStats)
#library(cowplot)
#library(hexbin)
#library("vsn")
#multivariate analysis (edgeR)
#if (!require("DT")) install.packages("DT"); library(DT) # for making interactive tables
#if (!require("plotly")) install.packages("plotly"); library(plotly) # for making interactive plots
#if (!require("gt")) install.packages("gt"); library(gt) # A layered 'grammar of tables' - think ggplot, but for tables
#if (!require("stargazer")) install.packages("stargazer"); library(stargazer)
#if (!require("rgl")) install.packages("rgl"); library(rgl)
#if (!require("dendextend")) install.packages("dendextend"); library(dendextend) # A way to customize dendrograms to make more complicated visualizations
#if (!require("ggpattern")) install.packages("ggpattern"); library(ggpattern)
#if (!require("data.table")) install.packages("data.table"); library(data.table)
#making pretty figs and tables
#devtools::install_github("thomasp85/patchwork")
#library(patchwork)
library(gt)
library("ggVennDiagram")
library(ggh4x)
library(patchwork)
library(ggpattern)
#lmerSeq analysis
#install dependencies
#install.packages(pkgs = c('pbapply', 'lmerTest', 'dplyr', 'devtools', 'magrittr', 'gtools', 'nlme', 'numDeriv', 'lavaSearch2'))
#download lmerSeq
#devtools::install_github("stop-pre16/lmerSeq", build_vignettes = T)
#other packages for lmerSeq
#if (!require("DESeq2")) install.packages("DESeq2");
#library(DESeq2)
#if (!require("magrittr")) install.packages("magrittr");
#library(magrittr)
#if (!require("limma")) install.packages("limma");
#library(limma)
#if (!require("dendextend")) install.packages("dendenxtend");
#library(dendextend) # To create dendrograms
#library(lmerSeq)
#for REBEL analysis
#install dependencies for RESCUE
#scater, scran, scuttle
#BiocManager::install("scater")
#BiocManager::install("scran")
#library(scran) #loads scuttle as part of it
#install Matrix.utils for REBEL
#remotes::install_github("cvarrichio/Matrix.utils")
#download RESCUE + REBEL - isntructions here: https://github.com/ewynn610/RESCUE and here: https://github.com/ewynn610/REBEL
#devtools::install_github("https://github.com/ewynn610/Rescue", build_vignettes = T)
#devtools::install_github("https://github.com/ewynn610/REBEL", build_vignettes = TRUE)
library(REBEL)
library(dplyr)
library(ggplot2)
library(plotly)
library(tidyr)
#phylogenies
library(ggtree)
library(ape)
library(treeio)
#general tools
library(stringr)Clear workspace
rm(list = ls())Read in data
Gene-level analysis
- Read in sample metadata
#read in sample metadata
meta <- read.delim("study_design.txt", header = TRUE, as.is = TRUE)
#make the SAMPLE column the rownames
rownames(meta) <- meta$SAMPLE
#add a names column (for lmerSeq later)
meta$names <- meta$SAMPLE
#capture sample labels from the design study file
sampleLabels <- meta$SAMPLE- Import variance-stabilized data
#import
vst_expr <- read_tsv("./variance_stabilized_counts_2025_03_12.tsv")
#convert to matrix
expr_mat <- as.matrix(vst_expr[,1:ncol(vst_expr)-1])
rownames(expr_mat) <- vst_expr$geneID- Import gene annotations
#get annotations for transcripts by loading file from ncbi ---- WAS EDITED TO MATCH ORs ADDED FROM Mitchell's ANNOTATION----
annotation_file <- ("Ppyr_annotations_ORmodified_2025_03_12.tsv")
#Tx <- read_tsv(annotation_file)
Tx <- read_tsv(annotation_file, col_names = c("target_id", "geneID", "name", "PpyrOR"), col_types = "cccc")
Tx <- as_tibble(Tx)Analysis: Which genes/ORs are most highly expressed?
- Calculate average expression across tissue * sex combinations
#calc average gene expression and standard deviation in relevant tissues, gene identity, OR identity
my.pheromone.data.df <- as_tibble(expr_mat) %>%
mutate(geneID = row.names(expr_mat)) %>%
rowwise() %>%
mutate(antenna.AVG = (SEL534ant + SEL535ant + SEL536ant + SEL543ant + SEL562ant + SEL627ant + SEL672ant)/7,
antenna.SD = sd(c(SEL534ant, SEL535ant, SEL536ant, SEL543ant, SEL562ant, SEL627ant, SEL672ant)),
leg.AVG = (SEL534BL + SEL535BL + SEL536BL + SEL543BL + SEL562BL + SEL627BL + SEL672BL)/7,
leg.SD = sd(c(SEL534BL, SEL535BL, SEL536BL, SEL543BL, SEL562BL, SEL627BL, SEL672BL)),
male.ant.AVG = (SEL534ant + SEL535ant + SEL536ant + SEL543ant)/4,
male.ant.SD = sd(c(SEL534ant, SEL535ant, SEL536ant, SEL543ant)),
male.leg.AVG = (SEL534BL + SEL535BL + SEL536BL + SEL543BL)/4,
male.leg.SD = sd(c(SEL534BL, SEL535BL, SEL536BL, SEL543BL)),
female.ant.AVG = (SEL562ant + SEL627ant + SEL672ant)/3,
female.ant.SD = sd(c(SEL562ant, SEL627ant, SEL672ant)),
female.leg.AVG = (SEL562BL + SEL627BL + SEL672BL)/3,
female.leg.SD = sd(c(SEL562BL, SEL627BL, SEL672BL)))
#ungroup the rows
my.pheromone.data.df <- my.pheromone.data.df %>%
ungroup()
#Add annotations (but only distinct geneIDs from the Tx file- has isoforms in it, and want to keep the PpyrOR name associated with the gene ID, even if it only matched 1 isoform)
Tx.with.PpyrORs.per.geneID <- Tx %>%
group_by(geneID) %>%
distinct(PpyrOR, .keep_all = TRUE) %>%
arrange(PpyrOR)
#checking
#nrow(Tx) #37237
#length(unique(Tx$geneID))#23,919
#nrow(Tx.with.PpyrORs.per.geneID) #23,934
#ungroup
Tx.with.PpyrORs.per.geneID.ungrouped <- Tx.with.PpyrORs.per.geneID %>%
ungroup()
#get one entry per geneID
Tx.one.entry.per.geneID<- distinct(Tx.with.PpyrORs.per.geneID.ungrouped, geneID, .keep_all = TRUE)
#checking
# nrow(Tx.one.entry.per.geneID) #23,919
# length(unique(Tx.one.entry.per.geneID$PpyrOR)) #103
my.pheromone.data.df.expanded <- left_join(my.pheromone.data.df, Tx.one.entry.per.geneID, by = "geneID")
#nrow(my.pheromone.data.df.expanded) #13,340
#add gene expression ranks by tissue/sex combo
my.pheromone.data.df.expanded$male.ant.RANK <- rank(-my.pheromone.data.df.expanded$male.ant.AVG, ties.method = "min")
my.pheromone.data.df.expanded$male.leg.RANK <- rank(-my.pheromone.data.df.expanded$male.leg.AVG, ties.method = "min")
my.pheromone.data.df.expanded$female.ant.RANK <- rank(-my.pheromone.data.df.expanded$female.ant.AVG, ties.method = "min")
my.pheromone.data.df.expanded$female.leg.RANK <- rank(-my.pheromone.data.df.expanded$female.leg.AVG, ties.method = "min")
#nrow(my.pheromone.data.df.expanded) 13,340
#get a subset of just ORs
my.OR.data.df.expanded <- my.pheromone.data.df.expanded %>%
filter(!is.na(PpyrOR))
#check to make sure is getting the "correct" subset
#nrow(my.OR.data.df.expanded) #45
### Extract the top 10 genes from each tissue * sex combo
#extract top 10
male.ant.top10 <- my.pheromone.data.df.expanded %>%
arrange(male.ant.RANK) %>%
select(name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)
#extract top 10
female.ant.top10 <-my.pheromone.data.df.expanded %>%
arrange(female.ant.RANK) %>%
select(name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)
#extract top 10
male.leg.top10 <- my.pheromone.data.df.expanded %>%
arrange(male.leg.RANK) %>%
select(geneID, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)
#extract top 10
female.leg.top10 <- my.pheromone.data.df.expanded %>%
arrange(female.leg.RANK) %>%
select(geneID, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)Which genes are most highly expressed?
- What are the top 10 genes in each tissue * sex condition?
Male antennae
#print nice table
male.ant.top10 |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
|||||
|---|---|---|---|---|---|---|---|---|
| Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | |
| PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA | 18.97 | 19.36 | 1 | 1 | 17.99 | 18.24 | 8 | 6 |
| PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA | 18.57 | 18.83 | 2 | 2 | 13.16 | 13.26 | 693 | 576 |
| PREDICTED: Photinus pyralis allergen Tha p 1-like (LOC116176659), mRNA | 17.99 | 18.05 | 3 | 3 | 17.82 | 17.66 | 10 | 12 |
| PREDICTED: Photinus pyralis uncharacterized LOC116179393 (LOC116179393), mRNA | 17.86 | 16.61 | 4 | 24 | 14.10 | 12.20 | 309 | 1,466 |
| PREDICTED: Photinus pyralis probable cytochrome P450 6a14 (LOC116162053), mRNA | 17.72 | 16.96 | 5 | 16 | 8.79 | 8.48 | 11,820 | 13,134 |
| PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA | 17.54 | 17.37 | 6 | 8 | 12.80 | 12.01 | 952 | 1,743 |
| PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA | 17.51 | 17.83 | 7 | 4 | 18.70 | 18.36 | 4 | 4 |
| PREDICTED: Photinus pyralis general odorant-binding protein 19d-like (LOC116167740), mRNA | 17.42 | 17.19 | 8 | 11 | 13.35 | 13.00 | 605 | 734 |
| PREDICTED: Photinus pyralis cytochrome P450 4C1-like (LOC116161249), mRNA | 17.20 | 17.67 | 9 | 6 | 15.66 | 15.62 | 75 | 68 |
| PREDICTED: Photinus pyralis uncharacterized LOC116181535 (LOC116181535), transcript variant X1, mRNA | 17.14 | 17.74 | 10 | 5 | 17.71 | 17.72 | 12 | 10 |
Female antennae
female.ant.top10 |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
|||||
|---|---|---|---|---|---|---|---|---|
| Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | |
| PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA | 18.97 | 19.36 | 1 | 1 | 17.99 | 18.24 | 8 | 6 |
| PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA | 18.57 | 18.83 | 2 | 2 | 13.16 | 13.26 | 693 | 576 |
| PREDICTED: Photinus pyralis allergen Tha p 1-like (LOC116176659), mRNA | 17.99 | 18.05 | 3 | 3 | 17.82 | 17.66 | 10 | 12 |
| PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA | 17.51 | 17.83 | 7 | 4 | 18.70 | 18.36 | 4 | 4 |
| PREDICTED: Photinus pyralis uncharacterized LOC116181535 (LOC116181535), transcript variant X1, mRNA | 17.14 | 17.74 | 10 | 5 | 17.71 | 17.72 | 12 | 10 |
| PREDICTED: Photinus pyralis cytochrome P450 4C1-like (LOC116161249), mRNA | 17.20 | 17.67 | 9 | 6 | 15.66 | 15.62 | 75 | 68 |
| PREDICTED: Photinus pyralis elongation factor 1-alpha (LOC116166585), transcript variant X1, mRNA | 16.69 | 17.44 | 17 | 7 | 17.16 | 17.50 | 19 | 14 |
| PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA | 17.54 | 17.37 | 6 | 8 | 12.80 | 12.01 | 952 | 1,743 |
| PREDICTED: Photinus pyralis uncharacterized LOC116174483 (LOC116174483), mRNA | 16.99 | 17.22 | 12 | 9 | 16.14 | 15.84 | 51 | 58 |
| PREDICTED: Photinus pyralis general odorant-binding protein 57d (LOC116162318), mRNA | 16.54 | 17.21 | 22 | 10 | 8.84 | 8.55 | 11,491 | 12,964 |
Male legs
male.leg.top10 |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| geneID | Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | ||
| 116178108 | PREDICTED: Photinus pyralis myosin heavy chain, muscle (LOC116178108), transcript variant X1, mRNA | 16.26 | 16.40 | 33 | 33 | 19.42 | 19.27 | 1 | 2 |
| 116180183 | PREDICTED: Photinus pyralis actin, muscle (LOC116180183), mRNA | 15.90 | 16.05 | 46 | 45 | 19.33 | 19.42 | 2 | 1 |
| 116168454 | PREDICTED: Photinus pyralis calcium-transporting ATPase sarcoplasmic/endoplasmic reticulum type (LOC116168454), transcript variant X1, mRNA | 16.51 | 16.37 | 23 | 35 | 19.05 | 18.55 | 3 | 3 |
| 116180227 | PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA | 17.51 | 17.83 | 7 | 4 | 18.70 | 18.36 | 4 | 4 |
| 116173867 | PREDICTED: Photinus pyralis paramyosin, long form-like (LOC116173867), mRNA | 15.04 | 14.92 | 90 | 110 | 18.52 | 18.19 | 5 | 7 |
| 116174773 | PREDICTED: Photinus pyralis titin (LOC116174773), mRNA | 14.41 | 14.53 | 166 | 176 | 18.40 | 18.25 | 6 | 5 |
| 116166781 | PREDICTED: Photinus pyralis troponin T, skeletal muscle (LOC116166781), transcript variant X1, mRNA | 14.93 | 15.14 | 99 | 91 | 18.01 | 17.91 | 7 | 8 |
| 116167561 | PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA | 18.97 | 19.36 | 1 | 1 | 17.99 | 18.24 | 8 | 6 |
| 116168198 | PREDICTED: Photinus pyralis muscle LIM protein Mlp84B-like (LOC116168198), transcript variant X1, mRNA | 13.83 | 13.65 | 304 | 399 | 17.84 | 17.66 | 9 | 11 |
| 116176659 | PREDICTED: Photinus pyralis allergen Tha p 1-like (LOC116176659), mRNA | 17.99 | 18.05 | 3 | 3 | 17.82 | 17.66 | 10 | 12 |
Female legs
female.leg.top10 |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| geneID | Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | ||
| 116180183 | PREDICTED: Photinus pyralis actin, muscle (LOC116180183), mRNA | 15.90 | 16.05 | 46 | 45 | 19.33 | 19.42 | 2 | 1 |
| 116178108 | PREDICTED: Photinus pyralis myosin heavy chain, muscle (LOC116178108), transcript variant X1, mRNA | 16.26 | 16.40 | 33 | 33 | 19.42 | 19.27 | 1 | 2 |
| 116168454 | PREDICTED: Photinus pyralis calcium-transporting ATPase sarcoplasmic/endoplasmic reticulum type (LOC116168454), transcript variant X1, mRNA | 16.51 | 16.37 | 23 | 35 | 19.05 | 18.55 | 3 | 3 |
| 116180227 | PREDICTED: Photinus pyralis ADP,ATP carrier protein (LOC116180227), mRNA | 17.51 | 17.83 | 7 | 4 | 18.70 | 18.36 | 4 | 4 |
| 116174773 | PREDICTED: Photinus pyralis titin (LOC116174773), mRNA | 14.41 | 14.53 | 166 | 176 | 18.40 | 18.25 | 6 | 5 |
| 116167561 | PREDICTED: Photinus pyralis cytochrome P450 4g15-like (LOC116167561), mRNA | 18.97 | 19.36 | 1 | 1 | 17.99 | 18.24 | 8 | 6 |
| 116173867 | PREDICTED: Photinus pyralis paramyosin, long form-like (LOC116173867), mRNA | 15.04 | 14.92 | 90 | 110 | 18.52 | 18.19 | 5 | 7 |
| 116166781 | PREDICTED: Photinus pyralis troponin T, skeletal muscle (LOC116166781), transcript variant X1, mRNA | 14.93 | 15.14 | 99 | 91 | 18.01 | 17.91 | 7 | 8 |
| 116175784 | PREDICTED: Photinus pyralis myosin regulatory light chain 2 (LOC116175784), mRNA | 14.71 | 14.78 | 118 | 132 | 17.78 | 17.73 | 11 | 9 |
| 116181535 | PREDICTED: Photinus pyralis uncharacterized LOC116181535 (LOC116181535), transcript variant X1, mRNA | 17.14 | 17.74 | 10 | 5 | 17.71 | 17.72 | 12 | 10 |
How many of the top ten genes are shared?
#display as a Venn Diagram
x = list(male.ant.top10$name, female.ant.top10$name, male.leg.top10$name, female.leg.top10$name)
ggVennDiagram(x,
category.names = c("Antenna, M", "Antenna, F", "Leg, M", "Leg, F")) + ggplot2::scale_fill_gradient(low = "#A0CBE8", high = "#4E79A7")
#sort(intersect(male.ant.top10$name, female.ant.top10$name))
#sort(intersect(male.leg.top10$name, female.leg.top10$name))
#sort(intersect(male.ant.top10$name, male.leg.top10$name))
#sort(intersect(female.ant.top10$name, female.leg.top10$name))- Overall, legs look similar to each other and different from antennae - most shared genes are actin, myosin, etc.
- Two genes are shared among the top 10 lists of all conditions - cytochrome P450 4g15-like (LOC116167561) and ADP,ATP carrier protein (LOC116180227).
Which ORs are most highly expressed?
What are the top 10 ORs in each tissue * sex condition?
- Orco not in top 10 for female legs, though present in all other lists (and first in both male and female antennae).
- PpyrOR18,23, 92 included in all top 10 lists.
- PpyrOR6, is highest male antenna-specific (in top 10 list).
- PpyrOR65 is antenna-specific (in top 10 lists).
- PpyrOR17PSE in top 10 female legs list (#10) - possibly functional?
### Extract the top 10 ORs from each tissue * sex combo
#extract top 10
male.ant.top10.OR <- my.OR.data.df.expanded %>%
arrange(male.ant.RANK) %>%
select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)
#extract top 10
female.ant.top10.OR <-my.OR.data.df.expanded %>%
arrange(female.ant.RANK) %>%
select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)
#extract top 10
male.leg.top10.OR <- my.OR.data.df.expanded %>%
arrange(male.leg.RANK) %>%
select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)
#extract top 10
female.leg.top10.OR <- my.OR.data.df.expanded %>%
arrange(female.leg.RANK) %>%
select(PpyrOR, name, male.ant.AVG, male.leg.AVG, female.ant.AVG, female.leg.AVG, male.ant.RANK, male.leg.RANK, female.ant.RANK, female.leg.RANK) %>%
slice_head(n = 10)Male antennae
#print nice table
male.ant.top10.OR |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| PpyrOR | Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | ||
| Ppyr/Orco | PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA | 13.49 | 13.09 | 429 | 670 | 8.68 | 8.41 | 12,528 | 13,230 |
| PpyrOR6 | PREDICTED: Photinus pyralis odorant receptor 94a-like (LOC116160183), transcript variant X1, mRNA | 11.84 | 9.14 | 1,954 | 10,155 | 8.36 | 8.37 | 13,291 | 13,257 |
| PpyrOR18 | PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA | 9.71 | 9.79 | 8,060 | 7,740 | 9.13 | 9.00 | 10,011 | 10,647 |
| PpyrOR23 | Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA | 9.69 | 10.06 | 8,149 | 6,818 | 9.00 | 8.98 | 10,609 | 10,747 |
| PpyrOR99 | PREDICTED: Photinus pyralis odorant receptor 43a-like (LOC116166584), mRNA | 9.69 | 9.14 | 8,165 | 10,154 | 8.49 | 8.47 | 13,171 | 13,153 |
| PpyrOR69 | PpyrOR69, NW_022170403 | 9.43 | 9.29 | 9,091 | 9,549 | 8.37 | 8.34 | 13,285 | 13,267 |
| PpyrOR92 | PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA | 9.42 | 9.84 | 9,108 | 7,579 | 8.71 | 8.65 | 12,316 | 12,593 |
| PpyrOR65 | PREDICTED: Photinus pyralis odorant receptor 67c-like (LOC116167462), mRNA | 9.36 | 9.57 | 9,318 | 8,499 | 8.41 | 8.52 | 13,259 | 13,045 |
| PpyrOR100 | PREDICTED: Photinus pyralis odorant receptor 43a-like (LOC116166593), mRNA | 9.31 | 8.77 | 9,544 | 12,091 | 8.44 | 8.50 | 13,238 | 13,097 |
| PpyrOR41 | PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116175472), mRNA | 9.30 | 9.33 | 9,578 | 9,426 | 8.37 | 8.38 | 13,288 | 13,249 |
Female antennae
#print nice table
female.ant.top10.OR |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| PpyrOR | Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | ||
| Ppyr/Orco | PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA | 13.49 | 13.09 | 429 | 670 | 8.68 | 8.41 | 12,528 | 13,230 |
| PpyrOR23 | Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA | 9.69 | 10.06 | 8,149 | 6,818 | 9.00 | 8.98 | 10,609 | 10,747 |
| PpyrOR92 | PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA | 9.42 | 9.84 | 9,108 | 7,579 | 8.71 | 8.65 | 12,316 | 12,593 |
| PpyrOR18 | PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA | 9.71 | 9.79 | 8,060 | 7,740 | 9.13 | 9.00 | 10,011 | 10,647 |
| PpyrOR28 | PREDICTED: Photinus pyralis odorant receptor 85c-like (LOC116176545), transcript variant X2, mRNA | 9.25 | 9.67 | 9,765 | 8,148 | 8.87 | 8.79 | 11,299 | 11,781 |
| PpyrOR88 | PpyrOR88, NW_022170878 | 9.19 | 9.59 | 10,004 | 8,445 | 8.34 | 8.34 | 13,297 | 13,267 |
| PpyrOR65 | PREDICTED: Photinus pyralis odorant receptor 67c-like (LOC116167462), mRNA | 9.36 | 9.57 | 9,318 | 8,499 | 8.41 | 8.52 | 13,259 | 13,045 |
| PpyrOR11 | PREDICTED: Photinus pyralis uncharacterized LOC116172672 (LOC116172672), mRNA | 9.11 | 9.53 | 10,271 | 8,672 | 8.34 | 8.34 | 13,297 | 13,267 |
| PpyrOR93 | PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161161), mRNA | 9.05 | 9.40 | 10,549 | 9,158 | 8.38 | 8.38 | 13,278 | 13,248 |
| PpyrOR70 | PpyrOR70, NW_022170403 | 9.19 | 9.38 | 9,983 | 9,232 | 8.44 | 8.34 | 13,239 | 13,267 |
Male legs
#print nice table
male.leg.top10.OR |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| PpyrOR | Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | ||
| PpyrOR18 | PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA | 9.71 | 9.79 | 8,060 | 7,740 | 9.13 | 9.00 | 10,011 | 10,647 |
| PpyrOR94 | PpyrOR94, NW_022170250 | 9.19 | 8.95 | 9,978 | 11,017 | 9.01 | 8.78 | 10,575 | 11,844 |
| PpyrOR23 | Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA | 9.69 | 10.06 | 8,149 | 6,818 | 9.00 | 8.98 | 10,609 | 10,747 |
| PpyrOR28 | PREDICTED: Photinus pyralis odorant receptor 85c-like (LOC116176545), transcript variant X2, mRNA | 9.25 | 9.67 | 9,765 | 8,148 | 8.87 | 8.79 | 11,299 | 11,781 |
| PpyrOR89 | PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116180111), transcript variant X1, mRNA | 8.84 | 9.10 | 11,687 | 10,292 | 8.84 | 9.00 | 11,483 | 10,654 |
| PpyrOR9 | Similar to XM_031483149.1 PREDICTED: Photinus pyralis odorant receptor Or2-like (LOC116167687), mRNA | 8.94 | 8.99 | 11,090 | 10,849 | 8.82 | 8.85 | 11,657 | 11,448 |
| PpyrOR25 | PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116177903), transcript variant X3, mRNA | 8.94 | 8.89 | 11,091 | 11,373 | 8.76 | 8.81 | 12,017 | 11,658 |
| PpyrOR92 | PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA | 9.42 | 9.84 | 9,108 | 7,579 | 8.71 | 8.65 | 12,316 | 12,593 |
| Ppyr/Orco | PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA | 13.49 | 13.09 | 429 | 670 | 8.68 | 8.41 | 12,528 | 13,230 |
| PpyrOR15 | PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116168902), mRNA | 9.00 | 9.05 | 10,795 | 10,545 | 8.67 | 8.38 | 12,587 | 13,251 |
Female legs
#print nice table
female.leg.top10.OR |>
gt() |>
tab_spanner(label = "Antenna", columns = contains("ant")) |>
tab_spanner(label = "Leg", columns = contains("leg")) |>
tab_spanner(label = "Mean Normalized Counts", columns = ends_with("AVG")) |>
tab_spanner(label = "Rank", columns = ends_with("RANK")) |>
cols_label(name = "Gene",
male.ant.AVG = "M",
male.leg.AVG = "M",
male.ant.RANK = "M",
male.leg.RANK = "M",
female.ant.AVG = "F",
female.leg.AVG = "F",
female.ant.RANK = "F",
female.leg.RANK = "F") |>
fmt_number(columns = ends_with("AVG"),
decimals = 2) |>
fmt_number(columns = ends_with("RANK"),
use_seps = TRUE,
decimals = 0) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = contains("leg") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#3C3C3C"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("male"))) |>
tab_style(style = list(cell_fill(color = "#717171"), cell_text(color = "white")),
locations = cells_body(columns = contains("ant") & starts_with("female"))) |>
tab_style(style = cell_fill(color = "#E5E5E5"),
locations = cells_body(columns = contains("leg") & starts_with("female"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Antenna", "Leg"))))|
Mean Normalized Counts
|
Rank
|
Mean Normalized Counts
|
Rank
|
||||||
|---|---|---|---|---|---|---|---|---|---|
| PpyrOR | Gene |
Antenna
|
Leg
|
||||||
| M | F | M | F | M | F | M | F | ||
| PpyrOR18 | PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178509), mRNA | 9.71 | 9.79 | 8,060 | 7,740 | 9.13 | 9.00 | 10,011 | 10,647 |
| PpyrOR89 | PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116180111), transcript variant X1, mRNA | 8.84 | 9.10 | 11,687 | 10,292 | 8.84 | 9.00 | 11,483 | 10,654 |
| PpyrOR23 | Similar to XM_031478590.1 PREDICTED: Photinus pyralis uncharacterized LOC116164413 (LOC116164413), transcript variant X1, mRNA | 9.69 | 10.06 | 8,149 | 6,818 | 9.00 | 8.98 | 10,609 | 10,747 |
| PpyrOR9 | Similar to XM_031483149.1 PREDICTED: Photinus pyralis odorant receptor Or2-like (LOC116167687), mRNA | 8.94 | 8.99 | 11,090 | 10,849 | 8.82 | 8.85 | 11,657 | 11,448 |
| PpyrOR25 | PREDICTED: Photinus pyralis odorant receptor 67a-like (LOC116177903), transcript variant X3, mRNA | 8.94 | 8.89 | 11,091 | 11,373 | 8.76 | 8.81 | 12,017 | 11,658 |
| PpyrOR28 | PREDICTED: Photinus pyralis odorant receptor 85c-like (LOC116176545), transcript variant X2, mRNA | 9.25 | 9.67 | 9,765 | 8,148 | 8.87 | 8.79 | 11,299 | 11,781 |
| PpyrOR94 | PpyrOR94, NW_022170250 | 9.19 | 8.95 | 9,978 | 11,017 | 9.01 | 8.78 | 10,575 | 11,844 |
| PpyrOR92 | PREDICTED: Photinus pyralis odorant receptor 82a-like (LOC116161158), transcript variant X2, mRNA | 9.42 | 9.84 | 9,108 | 7,579 | 8.71 | 8.65 | 12,316 | 12,593 |
| PpyrOR19 | PREDICTED: Photinus pyralis odorant receptor 4-like (LOC116178510), transcript variant X2, mRNA | 8.72 | 8.95 | 12,455 | 11,018 | 8.60 | 8.57 | 12,866 | 12,912 |
| PpyrOR17PSE | PpyrOR17PSE, NW_022170819 | 8.55 | 8.66 | 13,093 | 12,641 | 8.53 | 8.53 | 13,075 | 13,020 |
How many of the top ten ORs are shared?
#display as a Venn Diagram
top.10.OR.lists <- list(male.ant.top10.OR$PpyrOR, female.ant.top10.OR$PpyrOR, male.leg.top10.OR$PpyrOR, female.leg.top10.OR$PpyrOR)
ggVennDiagram(top.10.OR.lists,
category.names = c("Antenna, M", "Antenna, F", "Leg, M", "Leg, F")) + ggplot2::scale_fill_gradient(low = "#A0CBE8", high = "#4E79A7")
Overall, legs and antennae look different from each other, though sex difference is greater in antennae. If the ORs are shared between M + F antennae, then they are also usually shared with legs.
Which ORs are NOT expressed?
ORs.only <- Tx[!is.na(Tx$PpyrOR),]
#nrow(ORs.only) #103
ORs.not.expressed <- sort(setdiff(ORs.only$PpyrOR, my.OR.data.df.expanded$PpyrOR))
n_PSE <- length(grep("PSE", Tx$PpyrOR)) #15 PSEs
n_PSE.not.expressed <- length(grep("PSE", ORs.not.expressed)) #12 PSEs- 58 (56%) ORs are “not expressed” - they didn’t pass filtering in the expression dataset. They may be expressed in other life stages.
- 12 of the 15 (80%) pseudogenized ORs (PSE) are in the “not expressed” list.
- 3 (20%) OR PSEs did pass filtering and are expressed in adults, either in antennae or legs.
Analysis: Differential expression
- Fit the model
#set factors for analysis
meta$TISSUE = factor(meta$TISSUE, levels = c("leg", "antenna"))
meta$SEX = factor(meta$SEX, levels = c("F", "M"))
#convert to Summarized Experiment format
rse = SummarizedExperiment(assays = SimpleList(normcounts = as.matrix(expr_mat)), colData = meta)
#fit the model
res_reb = rebelFit(object = rse, fixedEffects = ~ TISSUE * SEX, #takes awhile
subjectVariable = "INDIVIDUAL",
pseudoBulk = T,
parallel = F,
outputFits = TRUE) - Calculate contrasts
#calculate DE for all contrasts
#male tissues: antenna vs leg
reb.mant_v_mleg <- rebelTest(res_reb, contrast = c(0, 1, 0, 1)) #male antennae vs leg
colnames(reb.mant_v_mleg) <- c("Estimate_mant_v_mleg", "Std.Error_mant_v_mleg", "t.value_mant_v_mleg", "df_mant_v_mleg", "p_val_raw_mant_v_mleg", "p_val_adj_mant_v_mleg")
reb.mant_v_mleg$geneID <- row.names(reb.mant_v_mleg) #add row names
#female tissues: antenna vs leg
reb.fant_v_fleg = rebelTest(res_reb, contrast = c(0, 1, 0, 0)) #female antennae vs leg - same as lmerseq?
colnames(reb.fant_v_fleg) <- c("Estimate_fant_v_fleg", "Std.Error_fant_v_fleg", "t.value_fant_v_fleg", "df_fant_v_fleg", "p_val_raw_fant_v_fleg", "p_val_adj_fant_v_fleg")
reb.fant_v_fleg$geneID <- row.names(reb.fant_v_fleg) #add row names
#antenna: male vs female
reb.mant_v_fant = rebelTest(res_reb, contrast = c(0, 0, 1, 1)) #male antennae vs female antenna - same as lmerseq?
colnames(reb.mant_v_fant) <- c("Estimate_mant_v_fant", "Std.Error_mant_v_fant", "t.value_mant_v_fant", "df_mant_v_fant", "p_val_raw_mant_v_fant", "p_val_adj_mant_v_fant")
reb.mant_v_fant$geneID <- row.names(reb.mant_v_fant) #add row names
#leg: male vs female
reb.mleg_v_fleg = rebelTest(res_reb, contrast = c(0, 0, 1, 0)) #male leg vs female leg - same as lmerseq?
colnames(reb.mleg_v_fleg) <- c("Estimate_mleg_v_fleg", "Std.Error_mleg_v_fleg", "t.value_mleg_v_fleg", "df_mleg_v_fleg", "p_val_raw_mleg_v_fleg", "p_val_adj_mleg_v_fleg")
reb.mleg_v_fleg$geneID <- row.names(reb.mleg_v_fleg) #add row names
#get into 1 big data table
#are all gene IDs in the same order? if so - can just cbind these together? #or just do join_by?
reb.tissue.comp <- left_join(reb.mant_v_mleg, reb.fant_v_fleg, by = "geneID")
reb.sex.comp <- left_join(reb.mant_v_fant, reb.mleg_v_fleg, by = "geneID")
reb.comp <- left_join(reb.tissue.comp, reb.sex.comp, by = "geneID")
#combine with my.pheromone.data.df.expanded (has annotations)
my.DE.data.df.expanded <- left_join(my.pheromone.data.df.expanded, reb.comp, by = "geneID")
#get ORs just on own
my.DE.OR.data.df.expanded <- my.DE.data.df.expanded %>%
filter(!is.na(PpyrOR))
#nrow(my.DE.OR.data.df.expanded) #45 - yes!
#export table (in nice order) as Supp table (big table)
#write.table(my.DE.data.df.expanded, "2025_03_12_DE_table_expanded.txt", row.names = FALSE, sep ='\t')- Collate and organize results
my.DE.data.df.expanded <- my.DE.data.df.expanded %>%
mutate(OR = case_when(!is.na(PpyrOR) ~ "OR",
TRUE ~"Other"),
DE_mant_v_mleg = ifelse(p_val_adj_mant_v_mleg < 0.05, "Significant", "Not significant"),
DE_fant_v_fleg = ifelse(p_val_adj_fant_v_fleg < 0.05, "Significant", "Not significant"),
DE_mant_v_fant = ifelse(p_val_adj_mant_v_fant < 0.05, "Significant", "Not significant"),
DE_mleg_v_fleg = ifelse(p_val_adj_mleg_v_fleg < 0.05, "Significant", "Not significant"),
Direction_mant_v_mleg = case_when(
# Estimate_mant_v_mleg > 0.6 & p_val_adj_mant_v_mleg < 0.05 ~ "Upregulated",
# Estimate_mant_v_mleg < -0.6 & p_val_adj_mant_v_mleg < 0.05 ~ "Downregulated",
Estimate_mant_v_mleg > 0 & p_val_adj_mant_v_mleg < 0.05 ~ "Upregulated",
Estimate_mant_v_mleg < 0 & p_val_adj_mant_v_mleg < 0.05 ~ "Downregulated",
TRUE ~ "Not significant"),
Direction_fant_v_fleg = case_when(
# Estimate_fant_v_fleg > 0.6 & p_val_adj_fant_v_fleg < 0.05 ~ "Upregulated",
# Estimate_fant_v_fleg < -0.6 & p_val_adj_fant_v_fleg < 0.05 ~ "Downregulated",
Estimate_fant_v_fleg > 0 & p_val_adj_fant_v_fleg < 0.05 ~ "Upregulated",
Estimate_fant_v_fleg < 0 & p_val_adj_fant_v_fleg < 0.05 ~ "Downregulated",
TRUE ~ "Not significant"),
Direction_mant_v_fant = case_when(
# Estimate_mant_v_fant > 0.6 & p_val_adj_mant_v_fant < 0.05 ~ "Upregulated",
# Estimate_mant_v_fant < -0.6 & p_val_adj_mant_v_fant < 0.05 ~ "Downregulated",
Estimate_mant_v_fant > 0 & p_val_adj_mant_v_fant < 0.05 ~ "Upregulated",
Estimate_mant_v_fant < 0 & p_val_adj_mant_v_fant < 0.05 ~ "Downregulated",
TRUE ~ "Not significant"),
Direction_mleg_v_fleg = case_when(
# Estimate_mleg_v_fleg > 0.6 & p_val_adj_mleg_v_fleg < 0.05 ~ "Upregulated",
# Estimate_mleg_v_fleg < -0.6 & p_val_adj_mleg_v_fleg < 0.05 ~ "Downregulated",
Estimate_mleg_v_fleg > 0 & p_val_adj_mleg_v_fleg < 0.05 ~ "Upregulated",
Estimate_mleg_v_fleg < 0 & p_val_adj_mleg_v_fleg < 0.05 ~ "Downregulated",
TRUE ~ "Not significant"))
my.DE.data.df.expanded$Direction_mant_v_mleg <- factor(my.DE.data.df.expanded$Direction_mant_v_mleg, levels = c("Upregulated", "Not significant", "Downregulated"))
my.DE.data.df.expanded$Direction_fant_v_fleg <- factor(my.DE.data.df.expanded$Direction_fant_v_fleg, levels = c("Upregulated", "Not significant", "Downregulated"))
my.DE.data.df.expanded$Direction_mant_v_fant <- factor(my.DE.data.df.expanded$Direction_mant_v_fant, levels = c("Upregulated", "Not significant", "Downregulated"))
my.DE.data.df.expanded$Direction_mleg_v_fleg <- factor(my.DE.data.df.expanded$Direction_mleg_v_fleg, levels = c("Upregulated", "Not significant", "Downregulated"))
#the OR table
my.DE.OR.data.df.expanded <- my.DE.data.df.expanded %>%
filter(!is.na(PpyrOR))
#Orco table
my.DE.ORco.data.df.expanded <- my.DE.data.df.expanded %>%
filter(PpyrOR == "Ppyr/Orco")Contrasts
Male: antenna v leg
How many genes and ORs are DE in male antenna vs leg?
m_tissue_DE <- my.DE.data.df.expanded %>%
filter(DE_mant_v_mleg == "Significant")
mant_v_mleg_DE_ORs <- my.DE.OR.data.df.expanded %>%
filter(DE_mant_v_mleg == "Significant") %>%
arrange(desc(Estimate_mant_v_mleg), p_val_adj_mant_v_mleg) %>%
select(PpyrOR, Direction_mant_v_mleg, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)- 3678 genes are significantly differentially expressed between male antenna and legs
- 32 ORs, out of 45 that passed the filter, (71.11)% are DE between male antenna and leg.
Which genes and ORs are they?
#for ggplot
my.mant_v_mleg.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), color = Direction_mant_v_mleg)) +
# geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) +
labs(color = "Expression",
x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated"))
#add circles, label Orco, OR6, male-specific ORs, and female-specific ORs
my.mant_v_mleg.plot.volcano.labels <- my.mant_v_mleg.plot.volcano +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_mleg, y=-log10(p_val_adj_mant_v_mleg), fill = Direction_mant_v_mleg), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#F28E2B", "#989ca3", "#4E79A7")) +
guides(fill = "none") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
#DE in male antennae
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
#male-specific DE
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), label = "OR100", color = "black", size = 3, hjust = -1.2, vjust = -1) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.15, yend = -log10(p_val_adj_mant_v_mleg) + .18), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), label = "OR67", color = "black", size = 3, hjust = -1.5, vjust = -0.9) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.7, yend = -log10(p_val_adj_mant_v_mleg) + 0.18), color = "grey28") +
#female-specific # F: 23, 28, 81, 15, 64, 59PSE
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), label = "OR59PSE", color = "black", size =3, hjust = -0.88, vjust = 1.1) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.8, yend = -log10(p_val_adj_mant_v_mleg) + -0.1), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), label = "OR23", color = "black", size = 3, hjust = -1.35, vjust = 1.37) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.4, yend = -log10(p_val_adj_mant_v_mleg) + -0.14), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), label = "OR28", color = "black", size =3, hjust = -1.52, vjust = 1.6) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.73, yend = -log10(p_val_adj_mant_v_mleg) + -0.17), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), label = "OR15", color = "black", size = 3, hjust = -1.56, vjust = 1.82) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.8, yend = -log10(p_val_adj_mant_v_mleg) + -0.2), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), label = "OR81", color = "black", size = 3, hjust = -1.65, vjust = 2.2) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.95, yend = -log10(p_val_adj_mant_v_mleg) + -0.25), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), label = "OR64", color = "black", size = 3, hjust = -1.64, vjust = 3) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), xend = Estimate_mant_v_mleg + 2.89, yend = -log10(p_val_adj_mant_v_mleg) + -0.35), color = "grey28") +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
my.mant_v_mleg.plot.volcano.labels
#for plotly
my.mant_v_mleg.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_mleg, y = -log10(p_val_adj_mant_v_mleg), color = Direction_mant_v_mleg, text = name)) +
geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) +
labs(color = "Male antenna v leg",
# x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?
my.mant_v_mleg.plot.circles <- my.mant_v_mleg.plot +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_mleg, y=-log10(p_val_adj_mant_v_mleg), fill = Direction_mant_v_mleg, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#F28E2B", "#989ca3", "#4E79A7")) +
guides(fill = "none")
ggplotly(my.mant_v_mleg.plot.circles)Circle with bold outline = ORs that passed filtering
Top 25 DE gene list
mant_v_mleg_DE_genes <- my.DE.data.df.expanded %>%
filter(p_val_adj_mant_v_mleg < 0.05) %>%
arrange(desc(Estimate_mant_v_mleg), p_val_adj_mant_v_mleg) %>%
select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
#nrow(mant_v_mleg_DE_genes) #3678 genes are DE
mant_v_mleg_DE_genes.top25 <- mant_v_mleg_DE_genes %>%
slice_head(n = 25)
mant_v_mleg_DE_genes.top25 |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(name = "Gene",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| Gene |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
| PREDICTED: Photinus pyralis probable cytochrome P450 6a14 (LOC116162053), mRNA | 8.93 | 0.00 | 8.48 | 0.00 | 0.77 | 0.23 | 0.31 | 0.46 |
| PREDICTED: Photinus pyralis general odorant-binding protein 57d (LOC116162318), mRNA | 7.69 | 0.00 | 8.66 | 0.00 | −0.68 | 0.13 | 0.29 | 0.37 |
| PREDICTED: Photinus pyralis protein takeout-like (LOC116181411), mRNA | 7.64 | 0.00 | 8.36 | 0.00 | −0.58 | 0.29 | 0.14 | 0.72 |
| PREDICTED: Photinus pyralis ejaculatory bulb-specific protein 3-like (LOC116176666), mRNA | 7.55 | 0.00 | 7.22 | 0.00 | 0.73 | 0.26 | 0.40 | 0.39 |
| PREDICTED: Photinus pyralis lipase 3-like (LOC116180894), transcript variant X1, mRNA | 7.48 | 0.00 | 7.45 | 0.00 | 0.12 | 0.82 | 0.09 | 0.73 |
| PREDICTED: Photinus pyralis uncharacterized LOC116165182 (LOC116165182), mRNA | 7.47 | 0.00 | 7.56 | 0.00 | 0.07 | 0.96 | 0.15 | 0.72 |
| PREDICTED: Photinus pyralis uncharacterized LOC116172495 (LOC116172495), mRNA | 7.37 | 0.00 | 7.14 | 0.00 | −0.36 | 0.51 | −0.59 | 0.16 |
| PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116170590), mRNA | 6.61 | 0.00 | 6.73 | 0.00 | 0.06 | 0.97 | 0.18 | 0.73 |
| PREDICTED: Photinus pyralis uncharacterized LOC116179743 (LOC116179743), transcript variant X1, mRNA | 6.28 | 0.00 | 4.95 | 0.00 | 1.74 | 0.11 | 0.41 | 0.54 |
| PREDICTED: Photinus pyralis general odorant-binding protein 72-like (LOC116170373), mRNA | 6.00 | 0.00 | 6.34 | 0.00 | −0.42 | 0.63 | −0.09 | 0.88 |
| PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA | 5.41 | 0.00 | 5.57 | 0.00 | −0.26 | 0.72 | −0.10 | 0.83 |
| PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161708), mRNA | 5.35 | 0.00 | 6.00 | 0.00 | −0.11 | 0.93 | 0.53 | 0.28 |
| PREDICTED: Photinus pyralis farnesol dehydrogenase-like (LOC116162028), mRNA | 5.32 | 0.00 | 5.17 | 0.00 | 0.25 | 0.75 | 0.10 | 0.83 |
| PREDICTED: Photinus pyralis probable cytochrome P450 6a13 (LOC116161618), mRNA | 5.00 | 0.00 | 4.91 | 0.00 | 0.12 | 0.96 | 0.02 | 0.97 |
| PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161490), transcript variant X1, mRNA | 4.93 | 0.00 | 4.38 | 0.00 | 1.12 | 0.50 | 0.57 | 0.55 |
| PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA | 4.82 | 0.00 | 4.68 | 0.00 | 0.40 | 0.51 | 0.27 | 0.46 |
| PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116164250), mRNA | 4.79 | 0.00 | 4.67 | 0.00 | 0.51 | 0.48 | 0.39 | 0.40 |
| PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA | 4.75 | 0.00 | 5.36 | 0.00 | 0.17 | 0.94 | 0.79 | 0.32 |
| PREDICTED: Photinus pyralis uncharacterized LOC116180967 (LOC116180967), mRNA | 4.30 | 0.00 | 4.35 | 0.00 | 0.22 | 0.90 | 0.27 | 0.69 |
| PREDICTED: Photinus pyralis odorant-binding protein 59a (LOC116179853), mRNA | 4.23 | 0.00 | 5.41 | 0.00 | 0.04 | 0.98 | 1.22 | 0.11 |
| PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161159), transcript variant X1, mRNA | 4.20 | 0.00 | 3.74 | 0.00 | 1.38 | 0.28 | 0.92 | 0.33 |
| PREDICTED: Photinus pyralis uncharacterized LOC116178615 (LOC116178615), mRNA | 4.20 | 0.00 | 4.67 | 0.00 | 0.21 | 0.86 | 0.69 | 0.28 |
| PREDICTED: Photinus pyralis glutathione S-transferase 1-like (LOC116167738), mRNA | 4.08 | 0.00 | 4.17 | 0.00 | 0.33 | 0.85 | 0.42 | 0.59 |
| PREDICTED: Photinus pyralis general odorant-binding protein 19d-like (LOC116167740), mRNA | 4.08 | 0.00 | 4.19 | 0.00 | 0.23 | 0.83 | 0.35 | 0.49 |
| PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116176983), transcript variant X1, mRNA | 3.97 | 0.00 | 4.20 | 0.00 | −0.04 | 0.98 | 0.19 | 0.73 |
DE OR list
mant_v_mleg_DE_ORs %>%
select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg) |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(PpyrOR = "OR",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| OR |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
| Ppyr/Orco | 4.82 | 0.00 | 4.68 | 0.00 | 0.40 | 0.51 | 0.27 | 0.46 |
| PpyrOR6 | 3.48 | 0.00 | 0.77 | 0.03 | 2.71 | 0.00 | −0.01 | 0.99 |
| PpyrOR99 | 1.20 | 0.00 | 0.67 | 0.01 | 0.55 | 0.13 | 0.02 | 0.95 |
| PpyrOR69 | 1.05 | 0.00 | 0.95 | 0.01 | 0.14 | 0.81 | 0.03 | 0.92 |
| PpyrOR65 | 0.95 | 0.01 | 1.05 | 0.01 | −0.21 | 0.69 | −0.11 | 0.73 |
| PpyrOR41 | 0.92 | 0.00 | 0.95 | 0.00 | −0.03 | 0.96 | −0.01 | 0.98 |
| PpyrOR100 | 0.87 | 0.00 | 0.27 | 0.07 | 0.54 | 0.09 | −0.06 | 0.76 |
| PpyrOR88 | 0.84 | 0.01 | 1.25 | 0.00 | −0.40 | 0.28 | 0.00 | 1.00 |
| PpyrOR11 | 0.77 | 0.01 | 1.19 | 0.01 | −0.42 | 0.30 | 0.00 | 1.00 |
| PpyrOR70 | 0.75 | 0.01 | 1.04 | 0.01 | −0.19 | 0.76 | 0.10 | 0.79 |
| PpyrOR38 | 0.73 | 0.00 | 0.67 | 0.00 | 0.06 | 0.90 | 0.00 | 1.00 |
| PpyrOR62 | 0.72 | 0.01 | 0.69 | 0.02 | 0.05 | 0.96 | 0.02 | 0.95 |
| PpyrOR92 | 0.71 | 0.02 | 1.19 | 0.01 | −0.42 | 0.41 | 0.06 | 0.87 |
| PpyrOR63 | 0.70 | 0.01 | 0.85 | 0.01 | −0.18 | 0.68 | −0.03 | 0.92 |
| PpyrOR93 | 0.67 | 0.00 | 1.02 | 0.00 | −0.35 | 0.19 | 0.00 | 1.00 |
| PpyrOR75 | 0.63 | 0.01 | 0.61 | 0.01 | 0.03 | 0.97 | 0.01 | 0.97 |
| PpyrOR24 | 0.60 | 0.02 | 0.99 | 0.01 | −0.42 | 0.32 | −0.03 | 0.94 |
| PpyrOR71 | 0.59 | 0.00 | 0.72 | 0.01 | −0.07 | 0.90 | 0.05 | 0.81 |
| PpyrOR68 | 0.59 | 0.01 | 0.48 | 0.02 | 0.20 | 0.58 | 0.09 | 0.69 |
| PpyrOR56PSE | 0.58 | 0.04 | 0.69 | 0.04 | −0.09 | 0.94 | 0.02 | 0.97 |
| PpyrOR18 | 0.58 | 0.01 | 0.80 | 0.01 | −0.08 | 0.92 | 0.13 | 0.66 |
| PpyrOR16 | 0.58 | 0.01 | 0.82 | 0.01 | −0.22 | 0.56 | 0.03 | 0.92 |
| PpyrOR98 | 0.57 | 0.02 | 0.85 | 0.01 | −0.26 | 0.68 | 0.02 | 0.97 |
| PpyrOR12 | 0.56 | 0.00 | 0.64 | 0.01 | −0.10 | 0.78 | −0.03 | 0.90 |
| PpyrOR76 | 0.55 | 0.01 | 0.77 | 0.00 | −0.30 | 0.42 | −0.08 | 0.75 |
| PpyrOR53 | 0.52 | 0.01 | 1.00 | 0.00 | −0.51 | 0.11 | −0.03 | 0.91 |
| PpyrOR58 | 0.49 | 0.02 | 0.59 | 0.02 | −0.07 | 0.93 | 0.03 | 0.91 |
| PpyrOR77 | 0.44 | 0.01 | 0.66 | 0.01 | −0.15 | 0.65 | 0.06 | 0.77 |
| PpyrOR61 | 0.42 | 0.01 | 0.38 | 0.03 | 0.04 | 0.96 | 0.00 | 1.00 |
| PpyrOR67 | 0.38 | 0.04 | 0.42 | 0.05 | −0.06 | 0.96 | −0.02 | 0.96 |
| PpyrOR66 | 0.35 | 0.02 | 0.51 | 0.01 | −0.24 | 0.36 | −0.08 | 0.64 |
| PpyrOR22 | 0.30 | 0.02 | 0.61 | 0.01 | −0.29 | 0.28 | 0.02 | 0.93 |
Female: antenna v leg
How many genes and ORs are DE in female antenna vs leg?
f_tissue_DE <- my.DE.data.df.expanded %>%
filter(DE_fant_v_fleg == "Significant")
fant_v_fleg_DE_ORs <- my.DE.OR.data.df.expanded %>%
filter(DE_fant_v_fleg == "Significant") %>%
arrange(desc(Estimate_fant_v_fleg), p_val_adj_fant_v_fleg) %>%
select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)- 3027 genes are significantly differentially expressed between female antenna and legs
- 36 ORs, out of 45 that passed the filter, (80)% are DE between female antenna and leg.
Which genes and ORs are they?
#for ggplot
my.fant_v_fleg.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), color = Direction_fant_v_fleg)) +
# geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) +
labs(color = "Female antenna v leg",
x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated"))
#add circles, label Orco, OR6, OR67 and 100, OR 15,23,28,59PSE,64,81
my.fant_v_fleg.plot.volcano.labels <- my.fant_v_fleg.plot.volcano +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_fant_v_fleg, y=-log10(p_val_adj_fant_v_fleg), fill = Direction_fant_v_fleg), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#F28E2B", "#989ca3", "#4E79A7")) +
guides(fill = "none") +
theme(legend.position = "none") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.05, vjust = -0.45) +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), label = "OR81", color = "black", size = 3, hjust = -2.22, vjust = -1.8) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR81"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4, yend = -log10(p_val_adj_fant_v_fleg) + 0.3), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), label = "OR28", color = "black", size = 3, hjust = -2.17, vjust = -3.6) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR28"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 3.9, yend = -log10(p_val_adj_fant_v_fleg) + 0.55), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), label = "OR59PSE", color = "black", size = 3, hjust = -1.35, vjust = -2.9) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR59PSE"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.26, yend = -log10(p_val_adj_fant_v_fleg) + 0.45), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -2.85, vjust = -2.4) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.05, yend = -log10(p_val_adj_fant_v_fleg) + 0.39), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), label = "OR64", color = "black", size = 3, hjust = -2.37, vjust = -1.8) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR64"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.25, yend = -log10(p_val_adj_fant_v_fleg) + 0.315), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), label = "OR23", color = "black", size = 3, hjust = -2.07, vjust = -0.6) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR23"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 3.7, yend = -log10(p_val_adj_fant_v_fleg) + 0.132), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), label = "OR15", color = "black", size = 3, hjust = -2.3, vjust = 0.4) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR15"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.15, yend = -log10(p_val_adj_fant_v_fleg) + 0.07), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), label = "OR67", color = "black", size = 3, hjust = -2.42, vjust = 1.4) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR67"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.36, yend = -log10(p_val_adj_fant_v_fleg) + -0.12), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), label = "OR100", color = "black", size = 3, hjust = -2.07, vjust = 1.78) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "PpyrOR100"), aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), xend = Estimate_fant_v_fleg + 4.5, yend = -log10(p_val_adj_fant_v_fleg) + -0.19), color = "grey28") +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
my.fant_v_fleg.plot.volcano.labels
#for plotly
my.fant_v_fleg.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_fant_v_fleg, y = -log10(p_val_adj_fant_v_fleg), color = Direction_fant_v_fleg, text = name)) +
geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) +
labs(color = "Female antenna v leg",
# x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?
my.fant_v_fleg.plot.circles <- my.fant_v_fleg.plot +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_fant_v_fleg, y=-log10(p_val_adj_fant_v_fleg), fill = Direction_fant_v_fleg, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#F28E2B", "#989ca3", "#4E79A7")) +
guides(fill = "none")
ggplotly(my.fant_v_fleg.plot.circles)Circle with bold outline = ORs that passed filtering
Top 25 DE gene list
fant_v_fleg_DE_genes <- my.DE.data.df.expanded %>%
filter(p_val_adj_fant_v_fleg < 0.05) %>%
arrange(desc(Estimate_fant_v_fleg), p_val_adj_fant_v_fleg) %>%
select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
fant_v_fleg_DE_genes.top25 <- fant_v_fleg_DE_genes %>%
slice_head(n = 25)
fant_v_fleg_DE_genes.top25 |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(name = "Gene",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| Gene |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
| PREDICTED: Photinus pyralis general odorant-binding protein 57d (LOC116162318), mRNA | 7.69 | 0.00 | 8.66 | 0.00 | −0.68 | 0.13 | 0.29 | 0.37 |
| PREDICTED: Photinus pyralis probable cytochrome P450 6a14 (LOC116162053), mRNA | 8.93 | 0.00 | 8.48 | 0.00 | 0.77 | 0.23 | 0.31 | 0.46 |
| PREDICTED: Photinus pyralis protein takeout-like (LOC116181411), mRNA | 7.64 | 0.00 | 8.36 | 0.00 | −0.58 | 0.29 | 0.14 | 0.72 |
| PREDICTED: Photinus pyralis uncharacterized LOC116165182 (LOC116165182), mRNA | 7.47 | 0.00 | 7.56 | 0.00 | 0.07 | 0.96 | 0.15 | 0.72 |
| PREDICTED: Photinus pyralis lipase 3-like (LOC116180894), transcript variant X1, mRNA | 7.48 | 0.00 | 7.45 | 0.00 | 0.12 | 0.82 | 0.09 | 0.73 |
| PREDICTED: Photinus pyralis ejaculatory bulb-specific protein 3-like (LOC116176666), mRNA | 7.55 | 0.00 | 7.22 | 0.00 | 0.73 | 0.26 | 0.40 | 0.39 |
| PREDICTED: Photinus pyralis uncharacterized LOC116172495 (LOC116172495), mRNA | 7.37 | 0.00 | 7.14 | 0.00 | −0.36 | 0.51 | −0.59 | 0.16 |
| PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116170590), mRNA | 6.61 | 0.00 | 6.73 | 0.00 | 0.06 | 0.97 | 0.18 | 0.73 |
| PREDICTED: Photinus pyralis general odorant-binding protein 72-like (LOC116170373), mRNA | 6.00 | 0.00 | 6.34 | 0.00 | −0.42 | 0.63 | −0.09 | 0.88 |
| PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161708), mRNA | 5.35 | 0.00 | 6.00 | 0.00 | −0.11 | 0.93 | 0.53 | 0.28 |
| PREDICTED: Photinus pyralis circadian clock-controlled protein-like (LOC116180925), mRNA | 5.41 | 0.00 | 5.57 | 0.00 | −0.26 | 0.72 | −0.10 | 0.83 |
| PREDICTED: Photinus pyralis odorant-binding protein 59a (LOC116179853), mRNA | 4.23 | 0.00 | 5.41 | 0.00 | 0.04 | 0.98 | 1.22 | 0.11 |
| PREDICTED: Photinus pyralis general odorant-binding protein 83a-like (LOC116171506), mRNA | 4.75 | 0.00 | 5.36 | 0.00 | 0.17 | 0.94 | 0.79 | 0.32 |
| PREDICTED: Photinus pyralis farnesol dehydrogenase-like (LOC116162028), mRNA | 5.32 | 0.00 | 5.17 | 0.00 | 0.25 | 0.75 | 0.10 | 0.83 |
| PREDICTED: Photinus pyralis uncharacterized LOC116179743 (LOC116179743), transcript variant X1, mRNA | 6.28 | 0.00 | 4.95 | 0.00 | 1.74 | 0.11 | 0.41 | 0.54 |
| PREDICTED: Photinus pyralis probable cytochrome P450 6a13 (LOC116161618), mRNA | 5.00 | 0.00 | 4.91 | 0.00 | 0.12 | 0.96 | 0.02 | 0.97 |
| PREDICTED: Photinus pyralis odorant receptor coreceptor (LOC116162435), mRNA | 4.82 | 0.00 | 4.68 | 0.00 | 0.40 | 0.51 | 0.27 | 0.46 |
| PREDICTED: Photinus pyralis uncharacterized LOC116178615 (LOC116178615), mRNA | 4.20 | 0.00 | 4.67 | 0.00 | 0.21 | 0.86 | 0.69 | 0.28 |
| PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116164250), mRNA | 4.79 | 0.00 | 4.67 | 0.00 | 0.51 | 0.48 | 0.39 | 0.40 |
| PREDICTED: Photinus pyralis uncharacterized LOC116172452 (LOC116172452), mRNA | 3.69 | 0.00 | 4.43 | 0.00 | −0.54 | 0.62 | 0.19 | 0.77 |
| PREDICTED: Photinus pyralis uncharacterized LOC116179393 (LOC116179393), mRNA | 3.76 | 0.01 | 4.41 | 0.01 | 1.25 | 0.47 | 1.90 | 0.17 |
| PREDICTED: Photinus pyralis alpha-tocopherol transfer protein-like (LOC116161490), transcript variant X1, mRNA | 4.93 | 0.00 | 4.38 | 0.00 | 1.12 | 0.50 | 0.57 | 0.55 |
| PREDICTED: Photinus pyralis uncharacterized LOC116180967 (LOC116180967), mRNA | 4.30 | 0.00 | 4.35 | 0.00 | 0.22 | 0.90 | 0.27 | 0.69 |
| PREDICTED: Photinus pyralis sensory neuron membrane protein 1-like (LOC116168667), mRNA | 3.94 | 0.00 | 4.26 | 0.00 | −0.14 | 0.92 | 0.18 | 0.71 |
| PREDICTED: Photinus pyralis medium-chain acyl-CoA ligase ACSF2, mitochondrial-like (LOC116176983), transcript variant X1, mRNA | 3.97 | 0.00 | 4.20 | 0.00 | −0.04 | 0.98 | 0.19 | 0.73 |
DE OR list
fant_v_fleg_DE_ORs |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(PpyrOR = "OR",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| OR |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
| Ppyr/Orco | 4.82 | 0.00 | 4.68 | 0.00 | 0.40 | 0.51 | 0.27 | 0.46 |
| PpyrOR88 | 0.84 | 0.01 | 1.25 | 0.00 | −0.40 | 0.28 | 0.00 | 1.00 |
| PpyrOR92 | 0.71 | 0.02 | 1.19 | 0.01 | −0.42 | 0.41 | 0.06 | 0.87 |
| PpyrOR11 | 0.77 | 0.01 | 1.19 | 0.01 | −0.42 | 0.30 | 0.00 | 1.00 |
| PpyrOR23 | 0.69 | 0.07 | 1.08 | 0.04 | −0.37 | 0.81 | 0.02 | 0.98 |
| PpyrOR65 | 0.95 | 0.01 | 1.05 | 0.01 | −0.21 | 0.69 | −0.11 | 0.73 |
| PpyrOR70 | 0.75 | 0.01 | 1.04 | 0.01 | −0.19 | 0.76 | 0.10 | 0.79 |
| PpyrOR93 | 0.67 | 0.00 | 1.02 | 0.00 | −0.35 | 0.19 | 0.00 | 1.00 |
| PpyrOR53 | 0.52 | 0.01 | 1.00 | 0.00 | −0.51 | 0.11 | −0.03 | 0.91 |
| PpyrOR24 | 0.60 | 0.02 | 0.99 | 0.01 | −0.42 | 0.32 | −0.03 | 0.94 |
| PpyrOR69 | 1.05 | 0.00 | 0.95 | 0.01 | 0.14 | 0.81 | 0.03 | 0.92 |
| PpyrOR41 | 0.92 | 0.00 | 0.95 | 0.00 | −0.03 | 0.96 | −0.01 | 0.98 |
| PpyrOR28 | 0.38 | 0.10 | 0.88 | 0.02 | −0.42 | 0.45 | 0.08 | 0.85 |
| PpyrOR63 | 0.70 | 0.01 | 0.85 | 0.01 | −0.18 | 0.68 | −0.03 | 0.92 |
| PpyrOR98 | 0.57 | 0.02 | 0.85 | 0.01 | −0.26 | 0.68 | 0.02 | 0.97 |
| PpyrOR16 | 0.58 | 0.01 | 0.82 | 0.01 | −0.22 | 0.56 | 0.03 | 0.92 |
| PpyrOR18 | 0.58 | 0.01 | 0.80 | 0.01 | −0.08 | 0.92 | 0.13 | 0.66 |
| PpyrOR76 | 0.55 | 0.01 | 0.77 | 0.00 | −0.30 | 0.42 | −0.08 | 0.75 |
| PpyrOR6 | 3.48 | 0.00 | 0.77 | 0.03 | 2.71 | 0.00 | −0.01 | 0.99 |
| PpyrOR81 | 0.17 | 0.18 | 0.77 | 0.01 | −0.63 | 0.09 | −0.03 | 0.90 |
| PpyrOR71 | 0.59 | 0.00 | 0.72 | 0.01 | −0.07 | 0.90 | 0.05 | 0.81 |
| PpyrOR56PSE | 0.58 | 0.04 | 0.69 | 0.04 | −0.09 | 0.94 | 0.02 | 0.97 |
| PpyrOR62 | 0.72 | 0.01 | 0.69 | 0.02 | 0.05 | 0.96 | 0.02 | 0.95 |
| PpyrOR38 | 0.73 | 0.00 | 0.67 | 0.00 | 0.06 | 0.90 | 0.00 | 1.00 |
| PpyrOR15 | 0.33 | 0.14 | 0.67 | 0.04 | −0.05 | 0.98 | 0.29 | 0.63 |
| PpyrOR99 | 1.20 | 0.00 | 0.67 | 0.01 | 0.55 | 0.13 | 0.02 | 0.95 |
| PpyrOR77 | 0.44 | 0.01 | 0.66 | 0.01 | −0.15 | 0.65 | 0.06 | 0.77 |
| PpyrOR12 | 0.56 | 0.00 | 0.64 | 0.01 | −0.10 | 0.78 | −0.03 | 0.90 |
| PpyrOR22 | 0.30 | 0.02 | 0.61 | 0.01 | −0.29 | 0.28 | 0.02 | 0.93 |
| PpyrOR75 | 0.63 | 0.01 | 0.61 | 0.01 | 0.03 | 0.97 | 0.01 | 0.97 |
| PpyrOR58 | 0.49 | 0.02 | 0.59 | 0.02 | −0.07 | 0.93 | 0.03 | 0.91 |
| PpyrOR64 | 0.21 | 0.21 | 0.53 | 0.04 | −0.36 | 0.55 | −0.04 | 0.94 |
| PpyrOR59PSE | 0.30 | 0.05 | 0.51 | 0.02 | −0.21 | 0.57 | 0.00 | 1.00 |
| PpyrOR66 | 0.35 | 0.02 | 0.51 | 0.01 | −0.24 | 0.36 | −0.08 | 0.64 |
| PpyrOR68 | 0.59 | 0.01 | 0.48 | 0.02 | 0.20 | 0.58 | 0.09 | 0.69 |
| PpyrOR61 | 0.42 | 0.01 | 0.38 | 0.03 | 0.04 | 0.96 | 0.00 | 1.00 |
Antennae: M v F
How many genes and ORs DE in male antenna vs female antenna?
antenna_DE <- my.DE.data.df.expanded %>%
filter(DE_mant_v_fant == "Significant")
mant_v_fant_DE_ORs <- my.DE.OR.data.df.expanded %>%
filter(DE_mant_v_fant == "Significant") %>%
arrange(desc(Estimate_mant_v_fant), p_val_adj_mant_v_fant) %>%
select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)- 2 genes are significantly differentially expressed between male and female antenna.
- 1 OR, out of 45 that passed the filter, (2.22)% is DE between male and female antenna.
Which genes and ORs are they?
#for ggplot
my.mant_v_fant.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_fant, y = -log10(p_val_adj_mant_v_fant), color = Direction_mant_v_fant)) +
# geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) +
labs(color = "Male v female antenna",
x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated"))
#add circles, label ORco and OR6
my.mant_v_fant.plot.volcano.labels <- my.mant_v_fant.plot.volcano +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_fant, y=-log10(p_val_adj_mant_v_fant), fill = Direction_mant_v_fant), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#F28E2B", "#989ca3", "#4E79A7")) +
guides(fill = "none") +
theme(legend.position = "none") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -1.7, vjust = 0.2) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), aes(x = Estimate_mant_v_fant, y = -log10(p_val_adj_mant_v_fant), xend = Estimate_mant_v_fant + 2.5, yend = -log10(p_val_adj_mant_v_fant) + 0.03), color = "grey28") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
my.mant_v_fant.plot.volcano.labels
#for plotly
my.mant_v_fant.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mant_v_fant, y = -log10(p_val_adj_mant_v_fant), color = Direction_mant_v_fant, text = name)) +
geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) +
labs(color = "Male v female antenna",
# x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?
my.mant_v_fant.plot.circles <- my.mant_v_fant.plot +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mant_v_fant, y=-log10(p_val_adj_mant_v_fant), fill = Direction_mant_v_fant, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#F28E2B", "#989ca3", "#4E79A7")) +
guides(fill = "none")
ggplotly(my.mant_v_fant.plot.circles)Circle with bold outline = ORs that passed filtering
DE gene list
mant_v_fant_DE_genes <- my.DE.data.df.expanded %>%
filter(p_val_adj_mant_v_fant < 0.05) %>%
arrange(desc(Estimate_mant_v_fant), p_val_adj_mant_v_fant) %>%
select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
mant_v_fant_DE_genes.top25 <- mant_v_fant_DE_genes %>%
slice_head(n = 25)
mant_v_fant_DE_genes.top25 |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(name = "Gene",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| Gene |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
| PREDICTED: Photinus pyralis odorant receptor 94a-like (LOC116160183), transcript variant X1, mRNA | 3.48 | 0.00 | 0.77 | 0.03 | 2.71 | 0.00 | −0.01 | 0.99 |
| PREDICTED: Photinus pyralis uncharacterized LOC116175393 (LOC116175393), transcript variant X1, mRNA | 1.21 | 0.01 | 0.00 | 1.00 | 1.38 | 0.05 | 0.16 | 0.67 |
DE OR list
mant_v_fant_DE_ORs <- my.DE.OR.data.df.expanded %>%
filter(DE_mant_v_fant == "Significant") %>%
arrange(desc(Estimate_mant_v_fant), p_val_adj_mant_v_fant) %>%
select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
mant_v_fant_DE_ORs |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(PpyrOR = "OR",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| OR |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
| PpyrOR6 | 3.48 | 0.00 | 0.77 | 0.03 | 2.71 | 0.00 | −0.01 | 0.99 |
Leg: M v F
How many genes DE in male leg vs female leg?
leg_DE <- my.DE.data.df.expanded %>%
filter(DE_mleg_v_fleg == "Significant")
mleg_v_fleg_DE_ORs <- my.DE.OR.data.df.expanded %>%
filter(DE_mleg_v_fleg == "Significant") %>%
arrange(desc(Estimate_mleg_v_fleg), p_val_adj_mleg_v_fleg) %>%
select(PpyrOR, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)- 10 genes are significantly differentially expressed between male and female legs
- 0 ORs, out of 45 that passed the filter, (0)% are DE between male and female legs.
Which genes are they?
#for ggplot
my.mleg_v_fleg.plot.volcano <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mleg_v_fleg, y = -log10(p_val_adj_mleg_v_fleg), color = Direction_mleg_v_fleg)) +
# geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-9, 9)) +
labs(color = "Male v female leg",
x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated"))
#add circles, label ORco and OR6
my.mleg_v_fleg.plot.volcano.labels <- my.mleg_v_fleg.plot.volcano +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mleg_v_fleg, y=-log10(p_val_adj_mleg_v_fleg), fill = Direction_mleg_v_fleg), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#989ca3", "#4E79A7", "#F28E2B")) +
guides(fill = "none") +
theme(legend.position = "none") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -2.2, vjust = 0.2) +
geom_segment(data = my.DE.OR.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), aes(x = Estimate_mleg_v_fleg, y = -log10(p_val_adj_mleg_v_fleg), xend = Estimate_mleg_v_fleg + 3.3, yend = -log10(p_val_adj_mleg_v_fleg) + 0.03), color = "grey28")+
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
my.mleg_v_fleg.plot.volcano.labels
#for plotly
my.mleg_v_fleg.plot <- ggplot(my.DE.data.df.expanded, aes(x = Estimate_mleg_v_fleg, y = -log10(p_val_adj_mleg_v_fleg), color = Direction_mleg_v_fleg, text = name)) +
geom_vline(xintercept = c(-0.6, 0.6), col = "gray", linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), col = "gray", linetype = "dashed") +
geom_point(size = 1) +
theme_set(theme_classic() +
theme(
axis.title.y = element_text(face = "bold", margin = margin(0,20,0,0), size = rel(1.1), color = 'black'),
axis.title.x = element_text(hjust = 0.5, face = "bold", margin = margin(20,0,0,0), size = rel(1.1), color = 'black'),
plot.title = element_text(hjust = 0.5)
)) +
coord_cartesian(ylim = c(0, 5), xlim = c(-10, 10)) +
labs(color = "Male v female leg",
# x = expression("log"[2]*"FC"), y = expression("-log"[10]*"p-value")) + # for ggplot
x = "log<sub>2</sub>FC", y = "-log<sub>10</sub>p-value") + #for plotly
scale_x_continuous(breaks = seq(-10, 10, 2)) +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"),
labels = c("Upregulated", "Not significant", "Downregulated")) #why doesn't this work?
my.mleg_v_fleg.plot.circles <- my.mleg_v_fleg.plot +
geom_point(data = my.DE.OR.data.df.expanded, mapping = aes(x=Estimate_mleg_v_fleg, y=-log10(p_val_adj_mleg_v_fleg), fill = Direction_mleg_v_fleg, text = PpyrOR), color = "black", size = 1.2, shape = 21) +
scale_fill_manual(values = c("#989ca3", "#4E79A7", "#F28E2B")) +
guides(fill = "none")
ggplotly(my.mleg_v_fleg.plot.circles)DE gene list
mleg_v_fleg_DE_genes <- my.DE.data.df.expanded %>%
filter(p_val_adj_mleg_v_fleg < 0.05) %>%
arrange(desc(Estimate_mleg_v_fleg), p_val_adj_mleg_v_fleg) %>%
select(name, Estimate_mant_v_mleg, p_val_adj_mant_v_mleg, Estimate_fant_v_fleg, p_val_adj_fant_v_fleg, Estimate_mant_v_fant, p_val_adj_mant_v_fant, Estimate_mleg_v_fleg, p_val_adj_mleg_v_fleg)
mleg_v_fleg_DE_genes.top25 <- mleg_v_fleg_DE_genes %>%
slice_head(n = 25)
mleg_v_fleg_DE_genes.top25 |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(name = "Gene",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| Gene |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
| PREDICTED: Photinus pyralis N-acylneuraminate-9-phosphatase (LOC116164423), mRNA | −1.36 | 0.00 | −0.20 | 0.23 | 0.25 | 0.63 | 1.40 | 0.03 |
| PREDICTED: Photinus pyralis uncharacterized LOC116178726 (LOC116178726), mRNA | −0.19 | 0.33 | 0.15 | 0.50 | 0.78 | 0.13 | 1.12 | 0.04 |
| PREDICTED: Photinus pyralis trifunctional nucleotide phosphoesterase protein YfkN (LOC116158805), transcript variant X1, mRNA | −0.43 | 0.04 | −0.18 | 0.35 | 0.76 | 0.10 | 1.01 | 0.04 |
| PREDICTED: Photinus pyralis uncharacterized LOC116174788 (LOC116174788), ncRNA | −1.02 | 0.01 | −0.10 | 0.60 | 0.06 | 0.95 | 0.98 | 0.04 |
| PREDICTED: Photinus pyralis mpv17-like protein (LOC116176952), mRNA | −0.33 | 0.07 | 0.00 | 1.00 | 0.51 | 0.20 | 0.84 | 0.04 |
| PREDICTED: Photinus pyralis titin homolog (LOC116173520), transcript variant X1, mRNA | 0.04 | 0.75 | 0.22 | 0.12 | 0.49 | 0.13 | 0.68 | 0.04 |
| PREDICTED: Photinus pyralis cysteine dioxygenase type 1 (LOC116159334), mRNA | 0.01 | 0.94 | −0.79 | 0.01 | 0.03 | 0.97 | −0.78 | 0.04 |
| PREDICTED: Photinus pyralis uncharacterized LOC116174750 (LOC116174750), transcript variant X1, mRNA | 0.02 | 0.87 | −0.05 | 0.64 | −1.42 | 0.09 | −1.49 | 0.04 |
| PREDICTED: Photinus pyralis phospholipase A1-like (LOC116181243), mRNA | 0.00 | 1.00 | −3.09 | 0.01 | −0.36 | 0.87 | −3.45 | 0.04 |
| PREDICTED: Photinus pyralis cytochrome P450 4c21-like (LOC116176575), mRNA | −0.31 | 0.16 | −0.99 | 0.01 | −2.78 | 0.09 | −3.47 | 0.04 |
DE OR list
mleg_v_fleg_DE_ORs |>
gt() |>
tab_spanner(label = "Male tissue AvL", columns = ends_with("mant_v_mleg")) |>
tab_spanner(label = "Female tissue AvL", columns = ends_with("fant_v_fleg")) |>
tab_spanner(label = "Antenna MvF", columns = ends_with("mant_v_fant")) |>
tab_spanner(label = "Leg MvF", columns = ends_with("mleg_v_fleg")) |>
cols_label(PpyrOR = "OR",
Estimate_mant_v_mleg = "log2FC",
p_val_adj_mant_v_mleg = "p",
Estimate_fant_v_fleg = "log2FC",
p_val_adj_fant_v_fleg = "p",
Estimate_mant_v_fant = "log2FC",
p_val_adj_mant_v_fant = "p",
Estimate_mleg_v_fleg = "log2FC",
p_val_adj_mleg_v_fleg = "p") |>
fmt_number(decimals = 2) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_mleg"))) |>
tab_style(style = cell_fill(color = "#D4D4D4"),
locations = cells_body(columns = ends_with("mant_v_fant"))) |>
tab_style(style = cell_text(weight = "bold"),
locations = list(cells_column_labels(), cells_column_spanners(spanners = c("Male tissue AvL", "Female tissue AvL", "Antenna MvF", "Leg MvF"))))| OR |
Male tissue AvL
|
Female tissue AvL
|
Antenna MvF
|
Leg MvF
|
||||
|---|---|---|---|---|---|---|---|---|
| log2FC | p | log2FC | p | log2FC | p | log2FC | p | |
Are upregulated ORs shared or unique?
x = list(mant_v_mleg_DE_ORs$PpyrOR, fant_v_fleg_DE_ORs$PpyrOR, mant_v_fant_DE_ORs$PpyrOR, mleg_v_fleg_DE_ORs$PpyrOR)
ggVennDiagram(x, category.names = c("Antenna v leg, M", "Antenna v leg, F", "Antenna, M v F", "Leg, M v F")) +
ggplot2::scale_fill_gradient(low = "#A0CBE8", high = "#4E79A7") + scale_x_continuous(expand = expansion(mult = .2))
Antenna-upregulated ORs unique to males
up_ant_in_males <- setdiff(mant_v_mleg_DE_ORs$PpyrOR, fant_v_fleg_DE_ORs$PpyrOR) #male unique
sort(up_ant_in_males)## [1] "PpyrOR100" "PpyrOR67"Antenna-upregulated ORs unique to females
up_ant_in_females <- setdiff(fant_v_fleg_DE_ORs$PpyrOR, mant_v_mleg_DE_ORs$PpyrOR) #female unique
sort(up_ant_in_females)## [1] "PpyrOR15" "PpyrOR23" "PpyrOR28" "PpyrOR59PSE" "PpyrOR64"
## [6] "PpyrOR81"Which ORs that are expressed aren’t DE?
not.in.males <- setdiff(my.DE.OR.data.df.expanded$PpyrOR, mant_v_mleg_DE_ORs$PpyrOR)
not.in.females <- setdiff(my.DE.OR.data.df.expanded$PpyrOR, fant_v_fleg_DE_ORs$PpyrOR)
try2<-my.DE.OR.data.df.expanded[which(my.DE.OR.data.df.expanded$PpyrOR %in% intersect(not.in.males, not.in.females)), ]
try2.long <- try2 %>% pivot_longer(cols = starts_with("SEL"),
names_to = "Sample",
values_to = "Expression")
try2.longer <- try2.long %>% mutate(Sex = case_when(Sample == "SEL534ant" ~ "M",
Sample == "SEL534BL" ~ "M",
Sample == "SEL535ant" ~ "M",
Sample == "SEL535BL" ~ "M",
Sample == "SEL536ant" ~ "M",
Sample == "SEL536BL" ~ "M",
Sample == "SEL543ant" ~ "M",
Sample == "SEL543BL" ~ "M",
Sample == "SEL562ant" ~ "F",
Sample == "SEL562BL" ~ "F",
Sample == "SEL627ant" ~ "F",
Sample == "SEL627BL" ~ "F",
Sample == "SEL672ant" ~ "F",
Sample == "SEL672BL" ~ "F",
TRUE ~ "not yet"))
try2.longer$Sex <- factor(try2.longer$Sex, levels = c("M", "F"))
#add tissue column
try2.longer <- try2.longer %>% mutate(Tissue = case_when(Sample == "SEL534ant" ~ "antenna",
Sample == "SEL534BL" ~ "leg",
Sample == "SEL535ant" ~ "antenna",
Sample == "SEL535BL" ~ "leg",
Sample == "SEL536ant" ~ "antenna",
Sample == "SEL536BL" ~ "leg",
Sample == "SEL543ant" ~ "antenna",
Sample == "SEL543BL" ~ "leg",
Sample == "SEL562ant" ~ "antenna",
Sample == "SEL562BL" ~ "leg",
Sample == "SEL627ant" ~ "antenna",
Sample == "SEL627BL" ~ "leg",
Sample == "SEL672ant" ~ "antenna",
Sample == "SEL672BL" ~ "leg",
TRUE ~ "not yet"))
#all.DE.ORs.idx.long$Tissue<- factor(all.DE.ORs.idx.long$Tissue, levels = c("antenna", "leg"))
#make sex*tissue column
try2.longer <- try2.longer %>% mutate(Treatment = paste(Sex, Tissue, sep = " "))
try2.longer$Treatment<- factor(try2.longer$Treatment, levels = c("M antenna", "F antenna", "M leg", "F leg"))
ggplot(try2.longer, aes(x = interaction(Sex, Tissue), y = Expression)) +
geom_boxplot(aes(color = Sex), outlier.shape = NA) + #outlier.size = 0.5
geom_jitter(size = 0.5, color = "black") +
facet_wrap2(~ PpyrOR, ncol = 2, scales = "free_y") + #strip = strip
guides(x = "axis_nested") +
# scale_fill_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
scale_color_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
# ylim(c(8,14.8)) +
ylab("Normalized counts") +
xlab("Sample") +
guides(color = "none")
- It seems like these genes have higher variance than others in the analysis.
- Conclusions:
- Many ORs are upregulated in antennae in both males (32) and females (36)
- PpyrOR6 is a top candidate for being involved in male-specific behavior (upregulated in antennae of both sexes, and upregulated in male antennae compared to females)
- PpyrOR100 and PpyrOR67 are also upregulated in male antennae vs. leg, but not females (and not sig. DE between sexes in antennae)
- PpyrOR15, PpyrOR23, PpyrOR28, PpyrOR59PSE, PpyrOR64, PpyrOR81 are upregulated in female antennae vs leg, but not males (and not sig. DE between sexes in antennae)
Testing model adequacy for DE ORs
# First, index the rows of the dataframe so that each gene has an index#
my.DE.data.df.expanded.idx <- my.DE.data.df.expanded %>%
mutate(index = 1:13340)
#nrow(my.DE.data.df.expanded.idx)
#then extract the DE ORs
all.DE.ORs.idx <- my.DE.data.df.expanded.idx %>%
filter(!is.na(PpyrOR)) %>%
filter(DE_mant_v_mleg == "Significant" | DE_fant_v_fleg == "Significant" | DE_mant_v_fant == "Significant" | DE_mleg_v_fleg == "Significant")
#nrow(all.DE.ORs.idx) #38 total DE
#then plot the fits
for(i in all.DE.ORs.idx$index){
print(plot(res_reb@fits[[i]], main = my.DE.data.df.expanded.idx$PpyrOR[i]))
}
Plots
ORs are generally upregulated in antennae
Distribution of DE estimates
my.DE.data.df.expanded.est.longer <- my.DE.data.df.expanded %>%
select(name, Estimate_mant_v_mleg, Estimate_fant_v_fleg, Estimate_mant_v_fant, Estimate_mleg_v_fleg, OR) %>%
pivot_longer(cols = starts_with("Estimate"), names_to = "Comparison", values_to = "Estimate")
my.DE.data.df.expanded.est.longer$Comparison <- factor(my.DE.data.df.expanded.est.longer$Comparison,
levels = c("Estimate_mant_v_mleg", "Estimate_fant_v_fleg", "Estimate_mant_v_fant", "Estimate_mleg_v_fleg"),
labels = c("M:AvL", "F:AvL", "A:MvF", "L:MvF"))
ggplot(my.DE.data.df.expanded.est.longer, aes(x=Estimate, fill = OR)) +
geom_histogram(alpha = 0.5, position = "identity") +
facet_grid(OR ~ Comparison, scales = "free_y") +
geom_vline(xintercept = 0, color = "gray", linetype = "dashed") +
guides(fill = "none")
- Note that there are no negative values for ORs in Antenna v Leg comparisons – the distribution is skewed right.
1:1 plots of leg vs antenna by sex
male.tissue.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=male.ant.AVG, x=male.leg.AVG, color = Direction_mant_v_mleg)) +
geom_point(shape=19, size=1, alpha = 0.75) +
ggtitle("Male") +
ylab("Antenna normalized counts") +
xlab("Leg normalized counts") +
theme_bw() +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Antenna-biased", "Not significant", "Leg-biased")) +
geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
theme_classic() +
theme(legend.title = element_blank())
male.tissue.1.2.1.labeled <- male.tissue.1.2.1 +
geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = male.ant.AVG, x = male.leg.AVG), shape = 21, color = "Black") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2)
female.tissue.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=female.ant.AVG, x=female.leg.AVG, color = Direction_fant_v_fleg)) +
geom_point(shape=19, size=1, alpha = 0.75) +
ggtitle("Female") +
ylab("Antenna normalized counts") +
xlab("Leg normalized counts") +
theme_bw() +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Antenna-biased", "Not significant", "Leg-biased")) +
geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
theme_classic() +
theme(legend.title = element_blank())
female.tissue.1.2.1.labeled <- female.tissue.1.2.1 +
geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = female.ant.AVG, x = female.leg.AVG), shape = 21, color = "Black") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2)
tissue.1.1.plots <- male.tissue.1.2.1.labeled / female.tissue.1.2.1.labeled +
plot_layout(guides = "collect", axis_titles = "collect") +
plot_annotation(tag_levels = 'A')
tissue.1.1.plots
#ggsave("2025_03_12_121_plot_tissue.png", plot = tissue.1.1.plots)PpyrOR6 is the only DE OR (one of two DE genes total) between male and female antennae
1:1 plot of male v female antennae (and legs)
ant.mant_v_fant.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=male.ant.AVG, x=female.ant.AVG, color = Direction_mant_v_fant)) +
geom_point(shape=19, size=1, alpha = 0.75) +
ggtitle("Antenna") +
ylab("Male normalized counts") +
xlab("Female normalized counts") +
theme_bw() +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Male-biased", "Not significant", "Female-biased")) +
geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
theme_classic() +
theme(legend.title = element_blank(), legend.position = "none")
ant.mant_v_fant.1.2.1.labeled <- ant.mant_v_fant.1.2.1 +
geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = male.ant.AVG, x = female.ant.AVG), shape = 21, color = "Black") +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "Ppyr/Orco"), label = "Orco", color = "black", size = 3, hjust = -0.2, vjust = 0.2) +
geom_text(data = my.DE.data.df.expanded %>% filter(PpyrOR == "PpyrOR6"), label = "OR6", color = "black", size = 3, hjust = -0.2, vjust = 0.2)
leg.mleg_v_fleg.1.2.1 <- ggplot(my.DE.data.df.expanded, aes(y=male.leg.AVG, x=female.leg.AVG, color = Direction_mleg_v_fleg)) +
geom_point(shape=19, size=1, alpha = 0.75) +
ggtitle("Leg") +
ylab("Male normalized counts") +
xlab("Female normalized counts") +
theme_bw() +
scale_color_manual(values = c("#F28E2B", "#989ca3", "#4E79A7"), labels = c("Male-biased", "Not significant", "Female-biased")) +
geom_abline(slope=1, intercept = 0, linetype = "dashed", color="grey28") +
theme_classic() +
theme(legend.title = element_blank())
leg.mleg_v_fleg.1.2.1.labeled <- leg.mleg_v_fleg.1.2.1 +
geom_point(data = my.DE.data.df.expanded %>% filter(!is.na(PpyrOR)), mapping = aes(y = male.leg.AVG, x = female.leg.AVG), shape = 21, color = "Black")
ant.mant_v_fant.1.2.1.labeled / leg.mleg_v_fleg.1.2.1.labeled +
plot_layout(guides = "collect", axis_titles = "collect") +
plot_annotation(tag_levels = 'A')
- The dotted line shows the 1:1 line
Volcano plots
volcanos.plot <- (my.mant_v_mleg.plot.volcano.labels + my.fant_v_fleg.plot.volcano.labels) / (my.mant_v_fant.plot.volcano.labels + my.mleg_v_fleg.plot.volcano.labels) +
plot_annotation(tag_levels = 'a') +
plot_layout(guides = "collect")
volcanos.plot
#ggsave("2025_03_12_volcanoes.svg", plot = volcanos.plot, device = "svg", width = 300, height = 300, dpi = 300, units = "mm")- Top left, male antennae vs leg
- Top right, female antennae vs leg
- Bottom left, male antennae vs female antennae
- Bottom right, male leg vs female leg
Plots of DE OR + Orco expression w/ p-values and sig bars/letters
Boxplot
all.DE.ORs.idx.long <- all.DE.ORs.idx %>% pivot_longer(cols = starts_with("SEL"),
names_to = "Sample",
values_to = "Expression")
#length(unique(all.DE.ORs.idx.long$PpyrOR))
#add sex column
all.DE.ORs.idx.long <- all.DE.ORs.idx.long %>% mutate(Sex = case_when(Sample == "SEL534ant" ~ "M",
Sample == "SEL534BL" ~ "M",
Sample == "SEL535ant" ~ "M",
Sample == "SEL535BL" ~ "M",
Sample == "SEL536ant" ~ "M",
Sample == "SEL536BL" ~ "M",
Sample == "SEL543ant" ~ "M",
Sample == "SEL543BL" ~ "M",
Sample == "SEL562ant" ~ "F",
Sample == "SEL562BL" ~ "F",
Sample == "SEL627ant" ~ "F",
Sample == "SEL627BL" ~ "F",
Sample == "SEL672ant" ~ "F",
Sample == "SEL672BL" ~ "F",
TRUE ~ "not yet"))
all.DE.ORs.idx.long$Sex <- factor(all.DE.ORs.idx.long$Sex, levels = c("M", "F"))
#add tissue column
all.DE.ORs.idx.long <- all.DE.ORs.idx.long %>% mutate(Tissue = case_when(Sample == "SEL534ant" ~ "antenna",
Sample == "SEL534BL" ~ "leg",
Sample == "SEL535ant" ~ "antenna",
Sample == "SEL535BL" ~ "leg",
Sample == "SEL536ant" ~ "antenna",
Sample == "SEL536BL" ~ "leg",
Sample == "SEL543ant" ~ "antenna",
Sample == "SEL543BL" ~ "leg",
Sample == "SEL562ant" ~ "antenna",
Sample == "SEL562BL" ~ "leg",
Sample == "SEL627ant" ~ "antenna",
Sample == "SEL627BL" ~ "leg",
Sample == "SEL672ant" ~ "antenna",
Sample == "SEL672BL" ~ "leg",
TRUE ~ "not yet"))
all.DE.ORs.idx.long$Tissue<- factor(all.DE.ORs.idx.long$Tissue, levels = c("antenna", "leg"))
#make sex*tissue column
all.DE.ORs.idx.long <- all.DE.ORs.idx.long %>% mutate(Treatment = paste(Sex, Tissue, sep = " "))
all.DE.ORs.idx.long$Treatment<- factor(all.DE.ORs.idx.long$Treatment, levels = c("M antenna", "F antenna", "M leg", "F leg"))
#make boxplots
df <- all.DE.ORs.idx.long
df$PpyrOR <- factor(df$PpyrOR, levels = c("Ppyr/Orco", "PpyrOR6", "PpyrOR11", "PpyrOR12", "PpyrOR15", "PpyrOR16", "PpyrOR18", "PpyrOR22", "PpyrOR23", "PpyrOR24", "PpyrOR28", "PpyrOR38", "PpyrOR41", "PpyrOR53", "PpyrOR56PSE", "PpyrOR58", "PpyrOR59PSE", "PpyrOR61", "PpyrOR62", "PpyrOR63", "PpyrOR64", "PpyrOR65", "PpyrOR66", "PpyrOR67", "PpyrOR68", "PpyrOR69", "PpyrOR70", "PpyrOR71", "PpyrOR75", "PpyrOR76", "PpyrOR77", "PpyrOR81", "PpyrOR88", "PpyrOR92", "PpyrOR93", "PpyrOR98", "PpyrOR99", "PpyrOR100"))
df$facet <- df$PpyrOR
rectangle_fill_sig_both <- "#bdbdbd"
rectangle_fill_sig_male_only <- "#636363"
rectangle_fill_sig_female_only <- "#f0f0f0"
strip <- strip_themed(background_x = elem_list_rect(fill = c(rectangle_fill_sig_both, #1, Orco *
rectangle_fill_sig_both, #2, Or6 *
rectangle_fill_sig_both, #3, Or11 *
rectangle_fill_sig_both, #4, Or12 *
rectangle_fill_sig_female_only, #5 , Or15 *
rectangle_fill_sig_both, #6, OR16 *
rectangle_fill_sig_both, #7, OR18 *
rectangle_fill_sig_both, #8, OR22 *
rectangle_fill_sig_female_only, #9, OR23 *
rectangle_fill_sig_both, #10, OR24 *
rectangle_fill_sig_female_only, #11, OR28 *
rectangle_fill_sig_both,#12, OR38 *
rectangle_fill_sig_both, # OR41 *
rectangle_fill_sig_both,#14, OR53 *
rectangle_fill_sig_both, # OR56PSE
rectangle_fill_sig_both,#16, OR58
rectangle_fill_sig_female_only, #, OR59PSE
rectangle_fill_sig_both,#18, OR61
rectangle_fill_sig_both, # 62
rectangle_fill_sig_both,#20, 63
rectangle_fill_sig_female_only, #, OR64
rectangle_fill_sig_both,#22, OR65
rectangle_fill_sig_both, #, OR66
rectangle_fill_sig_male_only,#24, OR67
rectangle_fill_sig_both, #, OR68
rectangle_fill_sig_both,#26, OR69
rectangle_fill_sig_both, # 70
rectangle_fill_sig_both,#28, 71
rectangle_fill_sig_both, # OR75
rectangle_fill_sig_both,#30, 76
rectangle_fill_sig_both, #, OR77
rectangle_fill_sig_female_only, #32, OR81
rectangle_fill_sig_both, #, OR88
rectangle_fill_sig_both,#34, 92
rectangle_fill_sig_both, #, 93
rectangle_fill_sig_both, #36, 98
rectangle_fill_sig_both, #99
rectangle_fill_sig_male_only), #38, 100
color = c(rep("black", 38)),
linewidth = c(1, 2, rep(1,36))),
text_x = elem_list_text(color = c(rep("black", 23), "white", rep("black", 13), "white")))#free scales
anno.free.Orco <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(14.2, 13.4),
y2 = c(14.5, 13.7),
xstar = c(2, 3),
ystar = c(14.6, 13.8),
lab = c("*", "*"),
PpyrOR = c("Ppyr/Orco", "Ppyr/Orco"))
anno.free.OR6 <- data.frame(x1 = c(1, 1, 2), #done
x2 = c(3, 2, 4),
y1 = c(13, 12.5, 12.0),
y2 = c(13.3, 12.8, 12.3),
xstar = c(2, 1.5, 3),
ystar = c(13.4, 12.9, 12.4),
lab = c("*", "*", "*"),
PpyrOR = c("PpyrOR6", "PpyrOR6", "PpyrOR6"))
anno.free.OR11 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(10.15, 9.9),
y2 = c(10.25, 10.0),
xstar = c(2, 3),
ystar = c(10.3, 10.05),
lab = c("*", "*"),
PpyrOR = c("PpyrOR11", "PpyrOR11"))
anno.free.OR12 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.26, 9.15),
y2 = c(9.3, 9.2),
xstar = c(2, 3),
ystar = c(9.31, 9.22),
lab = c("*", "*"),
PpyrOR = c("PpyrOR12", "PpyrOR12"))
anno.free.OR15 <- data.frame(x1 = c(2), #done, female only
x2 = c(4),
y1 = c(9.6),
y2 = c(9.65),
xstar = c(3),
ystar = c(9.67),
lab = c("*"),
PpyrOR = c("PpyrOR15"))
anno.free.OR16 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.72, 9.6),
y2 = c(9.77, 9.65),
xstar = c(2, 3),
ystar = c(9.79, 9.67),
lab = c("*", "*"),
PpyrOR = c("PpyrOR16", "PpyrOR16"))
anno.free.OR18 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(10.3, 10.1),
y2 = c(10.4, 10.2),
xstar = c(2, 3),
ystar = c(10.45, 10.22),
lab = c("*", "*"),
PpyrOR = c("PpyrOR18", "PpyrOR18"))
anno.free.OR22 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.3, 9.1),
y2 = c(9.4, 9.2),
xstar = c(2, 3),
ystar = c(9.45, 9.22),
lab = c("*", "*"),
PpyrOR = c("PpyrOR22", "PpyrOR22"))
anno.free.OR23 <- data.frame(x1 = c(2), #done, female only
x2 = c(4),
y1 = c(10.4),
y2 = c(10.5),
xstar = c(3),
ystar = c(10.6),
lab = c("*"),
PpyrOR = c("PpyrOR23"))
anno.free.OR24 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(10, 9.75),
y2 = c(10.1, 9.85),
xstar = c(2, 3),
ystar = c(10.15, 9.9),
lab = c("*", "*"),
PpyrOR = c("PpyrOR24", "PpyrOR24"))
anno.free.OR28 <- data.frame(x1 = c(2), #done, female only
x2 = c(4),
y1 = c(10.1),
y2 = c(10.18),
xstar = c(3),
ystar = c(10.2),
lab = c("*"),
PpyrOR = c("PpyrOR28"))
anno.free.OR38 <- data.frame(x1 = c(1, 2),
x2 = c(3, 4),
y1 = c(9.3, 9.1),
y2 = c(9.4, 9.2),
xstar = c(2, 3),
ystar = c(9.45, 9.22),
lab = c("*", "*"),
PpyrOR = c("PpyrOR38", "PpyrOR38"))
anno.free.OR41 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.6, 9.4),
y2 = c(9.7, 9.5),
xstar = c(2, 3),
ystar = c(9.75, 9.53),
lab = c("*", "*"),
PpyrOR = c("PpyrOR41", "PpyrOR41"))
anno.free.OR53 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.6, 9.4),
y2 = c(9.7, 9.5),
xstar = c(2, 3),
ystar = c(9.75, 9.53),
lab = c("*", "*"),
PpyrOR = c("PpyrOR53", "PpyrOR53"))
anno.free.OR56 <- data.frame(x1 = c(1, 2), #done, both
x2 = c(3, 4),
y1 = c(9.45, 9.2),
y2 = c(9.50, 9.25),
xstar = c(2, 3),
ystar = c(9.53, 9.28),
lab = c("*", "*"),
PpyrOR = c("PpyrOR56PSE", "PpyrOR56PSE"))
anno.free.OR58 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.4, 9.25),
y2 = c(9.45, 9.3),
xstar = c(2, 3),
ystar = c(9.48, 9.33),
lab = c("*", "*"),
PpyrOR = c("PpyrOR58", "PpyrOR58"))
anno.free.OR59 <- data.frame(x1 = c(2), #female only
x2 = c(4),
y1 = c(9),
y2 = c(9.05),
xstar = c(3),
ystar = c(9.065),
lab = c("*"),
PpyrOR = c("PpyrOR59PSE"))
anno.free.OR61 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9, 8.9),
y2 = c(9.05, 8.95),
xstar = c(2, 3),
ystar = c(9.065, 8.965),
lab = c("*", "*"),
PpyrOR = c("PpyrOR61", "PpyrOR61"))
anno.free.OR62 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.6, 9.4),
y2 = c(9.7, 9.5),
xstar = c(2, 3),
ystar = c(9.735, 9.535),
lab = c("*", "*"),
PpyrOR = c("PpyrOR62", "PpyrOR62"))
anno.free.OR63 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.6, 9.45),
y2 = c(9.65, 9.5),
xstar = c(2, 3),
ystar = c(9.7, 9.535),
lab = c("*", "*"),
PpyrOR = c("PpyrOR63", "PpyrOR63"))
anno.free.OR64 <- data.frame(x1 = c(2), #done, female only
x2 = c(4),
y1 = c(9.25),
y2 = c(9.3),
xstar = c(3),
ystar = c(9.335),
lab = c("*"),
PpyrOR = c("PpyrOR64"))
anno.free.OR65 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(10.0, 9.8),
y2 = c(10.1, 9.9),
xstar = c(2, 3),
ystar = c(10.135, 9.935),
lab = c("*", "*"),
PpyrOR = c("PpyrOR65", "PpyrOR65"))
anno.free.OR66 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.22, 9.1),
y2 = c(9.27, 9.15),
xstar = c(2, 3),
ystar = c(9.3, 9.18),
lab = c("*", "*"),
PpyrOR = c("PpyrOR66", "PpyrOR66"))
anno.free.OR67 <- data.frame(x1 = c(1), #done, male only
x2 = c(3),
y1 = c(9.1),
y2 = c(9.15),
xstar = c(2),
ystar = c(9.165),
lab = c("*", "*"),
PpyrOR = c("PpyrOR67"))
anno.free.OR68 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.3, 9.1),
y2 = c(9.4, 9.2),
xstar = c(2, 3),
ystar = c(9.45, 9.22),
lab = c("*", "*"),
PpyrOR = c("PpyrOR68", "PpyrOR68"))
anno.free.OR69 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.8, 9.55),
y2 = c(9.9, 9.65),
xstar = c(2, 3),
ystar = c(9.95, 9.7),
lab = c("*", "*"),
PpyrOR = c("PpyrOR69", "PpyrOR69"))
anno.free.OR70 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(10.1, 9.8),
y2 = c(10.2, 9.9),
xstar = c(2, 3),
ystar = c(10.25,9.95),
lab = c("*", "*"),
PpyrOR = c("PpyrOR70", "PpyrOR70"))
anno.free.OR71 <- data.frame(x1 = c(1, 2),
x2 = c(3, 4),
y1 = c(9.45, 9.3),
y2 = c(9.5, 9.35),
xstar = c(2, 3),
ystar = c(9.535, 9.385),
lab = c("*", "*"),
PpyrOR = c("PpyrOR71", "PpyrOR71"))
anno.free.OR75 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.55, 9.4),
y2 = c(9.6, 9.45),
xstar = c(2, 3),
ystar = c(9.635, 9.485),
lab = c("*", "*"),
PpyrOR = c("PpyrOR75", "PpyrOR75"))
anno.free.OR76 <- data.frame(x1 = c(1, 2),
x2 = c(3, 4),
y1 = c(9.6, 9.4),
y2 = c(9.65, 9.45),
xstar = c(2, 3),
ystar = c(9.7, 9.5),
lab = c("*", "*"),
PpyrOR = c("PpyrOR76", "PpyrOR76"))
anno.free.OR77 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.3, 9.1),
y2 = c(9.4, 9.2),
xstar = c(2, 3),
ystar = c(9.45, 9.22),
lab = c("*", "*"),
PpyrOR = c("PpyrOR77", "PpyrOR77"))
anno.free.OR81 <- data.frame(x1 = c(2), #done, female only
x2 = c(4),
y1 = c(9.5),
y2 = c(9.55),
xstar = c(3),
ystar = c(9.58),
lab = c("*"),
PpyrOR = c("PpyrOR81"))
anno.free.OR88 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(10.3, 10),
y2 = c(10.4, 10.1),
xstar = c(2, 3),
ystar = c(10.5, 10.2),
lab = c("*", "*"),
PpyrOR = c("PpyrOR88", "PpyrOR88"))
anno.free.OR92 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(10.6, 10.3),
y2 = c(10.7, 10.4),
xstar = c(2, 3),
ystar = c(10.8, 10.5),
lab = c("*", "*"),
PpyrOR = c("PpyrOR92", "PpyrOR92"))
anno.free.OR93 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.75, 9.6),
y2 = c(9.8, 9.65),
xstar = c(2, 3),
ystar = c(9.83, 9.68),
lab = c("*", "*"),
PpyrOR = c("PpyrOR93", "PpyrOR93"))
anno.free.OR98 <- data.frame(x1 = c(1, 2), #done
x2 = c(3, 4),
y1 = c(9.75, 9.6),
y2 = c(9.8, 9.65),
xstar = c(2, 3),
ystar = c(9.83, 9.68),
lab = c("*", "*"),
PpyrOR = c("PpyrOR98", "PpyrOR98"))
anno.free.OR99 <- data.frame(x1 = c(1, 2),
x2 = c(3, 4),
y1 = c(10.1, 9.8),
y2 = c(10.2, 9.9),
xstar = c(2, 3),
ystar = c(10.3, 10.0),
lab = c("*", "*"),
PpyrOR = c("PpyrOR99", "PpyrOR99"))
anno.free.OR100 <- data.frame(x1 = c(1), #done, male only
x2 = c(3),
y1 = c(9.45),
y2 = c(9.5),
xstar = c(2),
ystar = c(9.535),
lab = c("*", "*"),
PpyrOR = c("PpyrOR100"))
anno.free <- rbind(anno.free.Orco, anno.free.OR6, anno.free.OR11, anno.free.OR12, anno.free.OR15, anno.free.OR16, anno.free.OR18, anno.free.OR22, anno.free.OR23, anno.free.OR24, anno.free.OR28, anno.free.OR38, anno.free.OR41, anno.free.OR53, anno.free.OR56, anno.free.OR58, anno.free.OR59, anno.free.OR61, anno.free.OR62, anno.free.OR63, anno.free.OR64, anno.free.OR65, anno.free.OR66, anno.free.OR67, anno.free.OR68, anno.free.OR69, anno.free.OR70, anno.free.OR71, anno.free.OR75, anno.free.OR76, anno.free.OR77, anno.free.OR81, anno.free.OR88, anno.free.OR92, anno.free.OR93, anno.free.OR98, anno.free.OR99, anno.free.OR100)
anno.free$PpyrOR <- factor(anno.free$PpyrOR, levels = c("Ppyr/Orco", "PpyrOR6", "PpyrOR11", "PpyrOR12", "PpyrOR15", "PpyrOR16", "PpyrOR18", "PpyrOR22", "PpyrOR23", "PpyrOR24", "PpyrOR28", "PpyrOR38", "PpyrOR41", "PpyrOR53", "PpyrOR56PSE", "PpyrOR58", "PpyrOR59PSE", "PpyrOR61", "PpyrOR62", "PpyrOR63", "PpyrOR64", "PpyrOR65", "PpyrOR66", "PpyrOR67", "PpyrOR68", "PpyrOR69", "PpyrOR70", "PpyrOR71", "PpyrOR75", "PpyrOR76", "PpyrOR77", "PpyrOR81", "PpyrOR88", "PpyrOR92", "PpyrOR93", "PpyrOR98", "PpyrOR99", "PpyrOR100"))
boxplots.free <-ggplot(df, aes(x = interaction(Sex, Tissue), y = Expression)) +
geom_boxplot(aes(color = Sex), outlier.shape = NA) + #outlier.size = 0.5
geom_jitter(size = 0.5, color = "black") +
facet_wrap2(~ PpyrOR, ncol = 5, strip = strip, scales = "free_y") +
guides(x = "axis_nested") +
# scale_fill_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
scale_color_manual(values = c("#4E79A7", "#F28E2B")) + #colors: "#F28E2B", "#989ca3", "#4E79A7" ; grayscale:"#636363", "#f0f0f0"
# ylim(c(8,14.8)) +
ylab("Normalized counts") +
xlab("Sample") +
guides(color = "none") +
geom_text(data = anno.free, aes(x = xstar, y = ystar, label = lab)) +
geom_segment(data = anno.free, aes(x = x1, xend = x1,
y = y1, yend = y2),
colour = "black") +
geom_segment(data = anno.free, aes(x = x2, xend = x2,
y = y1, yend = y2),
colour = "black") +
geom_segment(data = anno.free, aes(x = x1, xend = x2,
y = y2, yend = y2),
colour = "black") +
xlab(element_blank())
boxplots.free
#ggsave("2025_03_12_boxplots_free.svg", plot = boxplots.free, device = "svg", width = 250, height = 350, dpi = 400, units = "mm")- * = p < 0.05
- blue = male, orange = female
- black jittered points = values for each sample
- shading of heading means significantly DE in male antennae vs legs (dark gray), female antennae v. legs (light gray), or both (medium gray)
- bold outline means significantly DE in male vs female antennae
Fold-change on linkage groups
#read in data file with this info
my.LG.data <- read.table("./bigger_OR_tibble_2025_03_12.txt", sep = "\t", header = TRUE)
#extract just one total length entry per OR
only.totals.tib <- my.LG.data %>% filter(Feature_title == "Total length")
#length(only.totals.tib$OR) #101
#make orco name compatible
#length(unique(my.DE.data.df.expanded$PpyrOR)) #46 - the ones that passed filtering
my.DE.data.df.expanded$PpyrOR <- gsub("Ppyr\\/Orco", "Ppyr\\\\Orco", my.DE.data.df.expanded$PpyrOR)
#join with expression dataframe
my.joined.tib <- left_join(only.totals.tib, my.DE.data.df.expanded, by = join_by(OR == PpyrOR))
#factor appropriately
my.joined.tib$OR_GROUP <- factor(my.joined.tib$OR_GROUP, levels = c("Orco", "OR1", "2A", "2B", "3", "4", "6"))
my.joined.tib$LG_common_name <- factor(my.joined.tib$LG_common_name, levels = c("LG1", "LG2", "LG3a (X)", "LG3b", "LG4", "LG5", "LG6", "LG7", "LG8", "LG9", "LG10"))
#my.joined.tib$suffix
my.joined.tib$suffix <- factor(my.joined.tib$suffix, levels = c("", "NTE", "PSE", "CTE"), labels = c("Full length", "Partial", "PSE", "Partial"))
#plot of LGs
#define option for non-scientific axis notation
options(scipen = 999)
#make tibble of LGs and format it (factor LGs)
LG_tibble <- tibble(LG_common_name = levels(my.joined.tib$LG_common_name), LG_length = c(70905644, 53657282, 22220000, 28664792, 50607672, 49173113, 47017841, 43653355, 37158542, 24427974, 20869423))
LG_tibble$LG_common_name <- factor(LG_tibble$LG_common_name, levels(my.joined.tib$LG_common_name))
my.joined.tib$DE_mant_v_mleg <- factor(my.joined.tib$DE_mant_v_mleg, levels = c("Significant", "Not significant", "NA"))
my.joined.tib$DE_fant_v_fleg <- factor(my.joined.tib$DE_fant_v_fleg, levels = c("Significant", "Not significant", "NA"))
my.joined.tib$DE_mant_v_fant <- factor(my.joined.tib$DE_mant_v_fant, levels = c("Significant", "Not significant", "NA"))
my.joined.tib$DE_mleg_v_fleg <- factor(my.joined.tib$DE_mleg_v_fleg, levels = c("Significant", "Not significant", "NA"))
my.joined.tib$Estimate_mant_v_mleg <- as.numeric(my.joined.tib$Estimate_mant_v_mleg)
my.joined.tib$Estimate_fant_v_fleg <- as.numeric(my.joined.tib$Estimate_fant_v_fleg)
my.joined.tib$Estimate_mant_v_fant <- as.numeric(my.joined.tib$Estimate_mant_v_fant)
my.joined.tib$Estimate_mleg_v_fleg <- as.numeric(my.joined.tib$Estimate_mleg_v_fleg)
library(scales)
#plot LG position
test <- my.joined.tib %>% select(OR, LG_common_name, converted_midpoint, Estimate_mant_v_mleg, Estimate_fant_v_fleg, Estimate_mant_v_fant, Estimate_mleg_v_fleg, DE_mant_v_mleg, DE_fant_v_fleg, DE_mant_v_fant, DE_mleg_v_fleg) %>%
pivot_longer(
cols = starts_with("Estimate") | starts_with("DE"),
cols_vary = "slowest",
names_to = c(".value", "Comparison"),
names_pattern = "(.+)_(...._v_....)")
#ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
# geom_bar(stat = "identity", color = "black", fill = "white") +
# theme_classic() +
# coord_flip() +
# scale_x_discrete(limits = rev) +
# scale_y_continuous(labels = scales::comma) +
# ylab("Length (bp)") +
# theme(axis.title.y = element_blank()) +
# geom_point(data = test, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate, color = DE, size = DE), shape = 21, alpha = 0.7, position = position_jitter(set.seed(42))) +
# scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
# scale_color_manual(values = c("black", "NA"), na.value = "gray48") +
# scale_size_manual(values = c(2, 2), na.value = 0.5) +
# facet_wrap(~ Comparison)
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_mleg, color = DE_mant_v_mleg, size = DE_mant_v_mleg), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
# geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_mleg, color = DE_mant_v_mleg, size = DE_mant_v_mleg), width = 0.4, height = 0, shape = 21, alpha = 0.7) +
scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
scale_color_manual(values = c("black", "gray48"), na.value = "gray70", name = "", labels = c("adj p \u003c 0.05", "adj p \u2265 0.05", "failed filtering")) +
scale_size_manual(values = c(2, 1), na.value = 0.5, name = "", labels = c("adj p \u003c 0.05", "adj p \u2265 0.05", "failed filtering")) +
# theme(legend.position = "none") +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Mant_v_mleg_fold-change_no_legend_on_LG_2025_03_12.svg", device = "svg")
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_fant_v_fleg, color = DE_fant_v_fleg, size = DE_fant_v_fleg), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
# geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_fant_v_fleg), size =2 , width = 0.4, height = 0, shape = 21, alpha = 0.6) +
scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC Female antennae v legs"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
theme(legend.position = "none") +
scale_color_manual(values = c("black", "gray48"), na.value = "gray70") +
scale_size_manual(values = c(2, 1), na.value = 0.5) +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Fant_v_fleg_fold-change_on_LG_2025_03_12.svg", device = "svg")
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_fant, color = DE_mant_v_fant, size = DE_mant_v_fant), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
# geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mant_v_fant), size =2 , width = 0.4, height = 0, shape = 21, alpha = 0.6) +
scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC Male v female antennae"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
theme(legend.position = "none") +
scale_color_manual(values = c("black", "gray48"), na.value = "gray70") +
scale_size_manual(values = c(2, 1), na.value = 0.5) +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Mant_v_fant_fold-change_on_LG_2025_03_12.svg", device = "svg")
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mleg_v_fleg, color = DE_mleg_v_fleg, size = DE_mleg_v_fleg), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
#geom_jitter(data = my.joined.tib, aes(x = LG_common_name, y = converted_midpoint, fill = Estimate_mleg_v_fleg), size =2 , width = 0.4, height = 0, shape = 21, alpha = 0.6) +
scale_fill_viridis_c(limits = c(0, 1.5), oob = squish, name = bquote(~log[2]~"FC Male v female legs"), na.value = "NA", labels = c("-0.7", "0.5", "1.0", "\u2265 1.5")) +
theme(legend.position = "none") +
scale_color_manual(values = c("black", "gray48"), na.value = "gray70") +
scale_size_manual(values = c(2, 1), na.value = 0.5) +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Mleg_v_fleg_fold-change_on_LG_2025_03_12.svg", device = "svg")
#for some reason, jittering changes when you assign the plot to a variable (like the set.seed doesn't work - so saved individually and combined in adobe)- Top to bottom:
- Male antennae v leg
- Female antennae v leg
- Male antennae v female antennae
- Male leg v female leg
Expression on linkage groups
#plot of LGs
#define option for non-scientific axis notation
options(scipen = 999)
#extract data for expression
expression_on_LGs.df <- my.joined.tib %>% select(OR, LG_common_name, converted_midpoint, male.ant.AVG, female.ant.AVG, male.leg.AVG, female.leg.AVG, DE_mant_v_mleg, DE_fant_v_fleg, DE_mant_v_fant, DE_mleg_v_fleg) %>%
pivot_longer(
cols = ends_with("AVG"),
cols_vary = "slowest", names_to = "Tissue", values_to = "Expression")
expression_on_LGs.df$Tissue <- factor(expression_on_LGs.df$Tissue, levels = c("male.ant.AVG", "female.ant.AVG", "male.leg.AVG", "female.leg.AVG"), labels = c("Male antennae", "Female antennae", "Male leg", "Female leg"))
expression_on_LGs.df$filtered <- "PASS"
expression_on_LGs.df$filtered[which(is.na(expression_on_LGs.df$Expression))] <- "FAIL"
expression_on_LGs.df$filtered <- factor(expression_on_LGs.df$filtered, levels = c("PASS", "FAIL"), labels = c("", "Filtered out"))
library(scales)
#plot LG position
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Male antennae") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Male_ant_norm_counts_on_LG_2025_03_12.svg", device = "svg")
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Female antennae") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
theme(legend.position = "none") +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Female_ant_norm_counts_on_LG_2025_03_12.svg", device = "svg")
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Male leg") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
theme(legend.position = "none") +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Male_leg_norm_counts_on_LG_2025_03_12.svg", device = "svg")
ggplot(LG_tibble, aes(x = LG_common_name, y = LG_length)) +
geom_bar(stat = "identity", color = "black", fill = "white") +
theme_classic() +
coord_flip() +
scale_x_discrete(limits = rev) +
scale_y_continuous(labels = scales::comma) +
ylab("Length (bp)") +
theme(axis.title.y = element_blank()) +
geom_point(data = expression_on_LGs.df %>% filter(Tissue == "Female leg") %>% select(OR, LG_common_name, converted_midpoint, Expression, filtered), aes(x = LG_common_name, y = converted_midpoint, fill = Expression, size = filtered), shape = 21, alpha = 0.7, position = position_jitter(set.seed(20))) +
scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", , labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
scale_size_manual(values = c(2, 0.5), name = "", labels = c("Passed filtering", "Failed filtering")) +
theme(legend.position = "none") +
force_panelsizes(rows = unit(85, "mm"),
cols = unit(85, "mm"))
#ggsave("Female_leg_norm_counts_on_LG_2025_03_12.svg", device = "svg")
#for some reason, jittering changes when you assign the plot to a variable (like the set.seed doesn't work - so saved individually and combined in adobe)- Top to bottom:
- Male antennae
- Female antennae
- Male leg
- Female leg
Expression on phylogeny
#load tree
#read in tree
tree <- read.tree("./run_2.contree_rooted_2025_03_11.newick")
#trim tree
#get the data from the tree
t_data <- as_tibble(tree)
#rename PpyrOR8
t_data <- t_data %>% mutate(label2 = label)
t_data$label2[which(t_data$label == "PpyrOR8b")] <- "PpyrOR8"
#rename Orco
t_data$label2[which(t_data$label == "Ppyr_Orco")] <- "PpyrOrco"
tree.rename <- rename_taxa(tree, t_data, label, label2)
t.rename_data <- as_tibble(tree.rename)
#just get the total OR info, 1 per OR, only firefly
only_totals <- my.joined.tib %>%
filter(feature == "total") %>%
filter(SPECIES == "Photinus pyralis")
#join the data
t_data <- left_join(t.rename_data, only_totals, by = join_by(label == OR_NAME))
#drop tips that are not Ppyr
to_drop <- which(is.na(t_data$LG))
tree_reduced <- drop.tip(tree.rename, to_drop)
#check tree looks good
#ggtree(tree_reduced) + geom_tiplab(size = 2) #good
#get reduced tree data
tr_data <- as_tibble(tree_reduced)
#add in OR info
tr_data <- left_join(tr_data, only_totals, by = join_by(label == OR_NAME))
tr_data$OR <- factor(tr_data$OR, levels = tr_data$OR)
tr_data$OR_GROUP <- factor(tr_data$OR_GROUP, levels = c("Orco", "OR1", "2A", "2B", "3", "4", "6"))
tr_data$LG_common_name <- factor(tr_data$LG_common_name, levels = c("LG1", "LG2", "LG3a (X)", "LG3b", "LG4", "LG5", "LG6", "LG7", "LG8", "LG9", "LG10"))
expression.tibble <- tr_data %>% select(parent, node, branch.length, label, male.ant.AVG, female.ant.AVG, male.leg.AVG, female.leg.AVG)
expression.tibble.longer <- expression.tibble %>% pivot_longer(ends_with("AVG"), names_to = "Sex x Tissue", values_to = "Normalized counts")
expression.tibble.longer <- expression.tibble.longer %>% mutate(Sex = ifelse(str_detect(`Sex x Tissue`, "female"), "Female", "Male"),
Tissue = ifelse(str_detect(`Sex x Tissue`, "ant"), "Antennae", "Legs"))
expression.tibble.longer$Sex <- factor(expression.tibble.longer$Sex, levels = c("Male", "Female"), labels = c("M", "F"))
expression.tibble.longer$Tissue <- factor(expression.tibble.longer$Tissue, levels = c("Antennae", "Legs"))
#orient tree
#attach expression data
#plot - can I have a series of circles @ end with expression values instead of 4 different trees? Or a heatmap?
library(aplot)
p <- ggtree(tree_reduced) %<+% tr_data
#add tip labels, group designations (color), and node labels to see which nodes to flip
#p + geom_tiplab(size = 2) +
# geom_tippoint(aes(fill = OR_GROUP), color = "black", shape = 21, size = 2) +
# geom_text(aes(label = node), size = 2, vjust = -1, hjust = 1.2)
# scale_fill_manual(values = c("#4E79A7","#E15759", "#B6992D", "#F28E2B", "#B07AA1", "#59A14F", "#797067"))
#flip the 2B node to be next to 2A for ease of visualization
p2 <- flip(p, 124, 107)
#plot flipped tree with labels and colors to check
#p2 +
# geom_tiplab(size = 2, hjust = -0.1) +
# geom_text(aes(label = node), size = 2, vjust = -1, hjust = 1.2) +
# geom_tippoint(aes(color = OR_GROUP), size = 1.5, alpha = 0.9)# +
# scale_color_manual(values = c("#4E79A7","#E15759", "#B6992D", "#F28E2B", "#B07AA1", "#59A14F", "#797067")) +
# theme(legend.position = "right") +
# xlim(0,10) +
# ylim(0, 101)
#get node values for labels
Group2A.MRCA <- getMRCA(tree_reduced, c("PpyrOR6", "PpyrOR2")) #group2A
Group2B.MRCA <- getMRCA(tree_reduced, c("PpyrOR81", "PpyrOR65")) #group2B
Group3.MRCA <-getMRCA(tree_reduced, c("PpyrOR92", "PpyrOR90")) #group3
Group4.MRCA <- getMRCA(tree_reduced, c("PpyrOR84", "PpyrOR87")) #group4
Group6.MRCA <- getMRCA(tree_reduced, c("PpyrOR100", "PpyrOR99")) #group6
z<-p2 +
geom_tiplab(align=TRUE, color = "white") + hexpand(0.15) + xlim(0,10) +
geom_cladelab(node=Group2A.MRCA, label="2A", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
geom_cladelab(node=Group2B.MRCA, label="2B", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
geom_cladelab(node=Group3.MRCA, label="3", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
geom_cladelab(node=Group4.MRCA, label="4", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
geom_cladelab(node=Group6.MRCA, label="6", align = TRUE, offset=0.9, barsize=1, fontsize=7) +
geom_strip("PpyrOR1", "PpyrOR1", label="OR1", align = TRUE, offset=-1.0, barsize=1, fontsize=3) +
annotate("segment", x = 3.45, y = 8.7, xend = 4.2, yend = 8.7, color = "gray48", lty = "1111", linewidth = 0.5)
z2 <- ggplot(expression.tibble.longer, aes(x = interaction(Sex, Tissue), y = label)) +
geom_tile(aes(fill = `Normalized counts`), na.rm = TRUE) +
scale_fill_viridis_c(limits = c(8.3, 11), oob = squish, name = bquote("Normalized counts"), na.value = "NA", labels = c("8.5", "9.0", "9.5", "10.0", "10.5", "\u2265 11.0")) +
theme_classic() +
xlab(NULL) +
ylab(NULL) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
guides(x = "axis_nested")
#+
# force_panelsizes(rows = unit(85, "mm"),
# cols = unit(45, "mm"))
z2 %>% insert_left(z)
#ggsave("Expression_on_phylo_2025_03_12.svg", device = "svg")Session information
sessionInfo()## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.7.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] aplot_0.2.3 scales_1.3.0
## [3] stringr_1.5.1 treeio_1.28.0
## [5] ape_5.8 ggtree_3.12.0
## [7] tidyr_1.3.1 plotly_4.10.4
## [9] dplyr_1.1.4 REBEL_0.1.0
## [11] ggpattern_1.1.1 patchwork_1.3.0.9000
## [13] ggh4x_0.2.8 ggplot2_3.5.1
## [15] ggVennDiagram_1.5.2 gt_0.11.1
## [17] readr_2.1.5 SummarizedExperiment_1.34.0
## [19] Biobase_2.64.0 GenomicRanges_1.56.2
## [21] GenomeInfoDb_1.40.1 IRanges_2.38.1
## [23] S4Vectors_0.42.1 BiocGenerics_0.50.0
## [25] MatrixGenerics_1.16.0 matrixStats_1.4.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2 farver_2.1.2
## [4] fastmap_1.2.0 lazyeval_0.2.2 digest_0.6.37
## [7] lifecycle_1.0.4 tidytree_0.4.6 magrittr_2.0.3
## [10] compiler_4.4.1 rlang_1.1.4 sass_0.4.9
## [13] tools_4.4.1 yaml_2.3.10 data.table_1.16.4
## [16] knitr_1.49 labeling_0.4.3 S4Arrays_1.4.1
## [19] htmlwidgets_1.6.4 bit_4.5.0.1 DelayedArray_0.30.1
## [22] xml2_1.3.6 abind_1.4-8 withr_3.0.2
## [25] purrr_1.0.2 grid_4.4.1 colorspace_2.1-1
## [28] cli_3.6.3 rmarkdown_2.29 crayon_1.5.3
## [31] generics_0.1.3 rstudioapi_0.17.1 httr_1.4.7
## [34] tzdb_0.4.0 cachem_1.1.0 zlibbioc_1.50.0
## [37] parallel_4.4.1 ggplotify_0.1.2 XVector_0.44.0
## [40] vctrs_0.6.5 yulab.utils_0.1.8 Matrix_1.7-0
## [43] jsonlite_1.8.9 gridGraphics_0.5-1 hms_1.1.3
## [46] bit64_4.5.2 crosstalk_1.2.1 jquerylib_0.1.4
## [49] glue_1.8.0 stringi_1.8.4 gtable_0.3.5
## [52] UCSC.utils_1.0.0 munsell_0.5.1 tibble_3.2.1
## [55] pillar_1.10.1 htmltools_0.5.8.1 GenomeInfoDbData_1.2.12
## [58] R6_2.5.1 vroom_1.6.5 evaluate_1.0.3
## [61] lattice_0.22-6 ggfun_0.1.8 bslib_0.8.0
## [64] Rcpp_1.0.14 SparseArray_1.4.8 nlme_3.1-166
## [67] xfun_0.50 fs_1.6.5 pkgconfig_2.0.3